X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=trimoligos.cpp;h=8f4cbe98a7488b9d552f07b6ffd17aeb5e858c23;hb=e10c72304ee071c0c40e0218a06d89dc4731cbc2;hp=cf818707c672c82acfdbd3bf896c9c0c2edbae36;hpb=d9c3ceb4147b827719046c7d742df14288a07722;p=mothur.git diff --git a/trimoligos.cpp b/trimoligos.cpp index cf81870..8f4cbe9 100644 --- a/trimoligos.cpp +++ b/trimoligos.cpp @@ -36,6 +36,27 @@ TrimOligos::TrimOligos(int p, int b, int l, int s, map pr, map pr, map br, vector lk, vector sp){ + try { + m = MothurOut::getInstance(); + + pdiffs = p; + bdiffs = b; + ldiffs = l; + sdiffs = s; + + ibarcodes = br; + iprimers = pr; + linker = lk; + spacer = sp; + } + catch(exception& e) { + m->errorOut(e, "TrimOligos", "TrimOligos"); + exit(1); + } +} +/********************************************************************/ +//strip, pdiffs, bdiffs, primers, barcodes, revPrimers TrimOligos::TrimOligos(int p, int b, map pr, map br, vector r){ try { m = MothurOut::getInstance(); @@ -129,8 +150,7 @@ int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){ } oligo = oligo.substr(0,alnLength); temp = temp.substr(0,alnLength); - - int numDiff = countDiffs(oligo, temp); + int numDiff = countDiffs(oligo, temp); if(numDiff < minDiff){ minDiff = numDiff; @@ -172,6 +192,177 @@ int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){ m->errorOut(e, "TrimOligos", "stripBarcode"); exit(1); } +} +//*******************************************************************/ +int TrimOligos::stripBarcode(Sequence& forwardSeq, Sequence& reverseSeq, QualityScores& forwardQual, QualityScores& reverseQual, int& group){ + try { + //look for forward barcode + string rawFSequence = forwardSeq.getUnaligned(); + string rawRSequence = reverseSeq.getUnaligned(); + int success = bdiffs + 1; //guilty until proven innocent + + //can you find the forward barcode + for(map::iterator it=ibarcodes.begin();it!=ibarcodes.end();it++){ + string foligo = it->second.forward; + string roligo = it->second.reverse; + + if(rawFSequence.length() < foligo.length()){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out + break; + } + if(rawRSequence.length() < roligo.length()){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out + break; + } + + if((compareDNASeq(foligo, rawFSequence.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawRSequence.substr((rawRSequence.length()-roligo.length()),roligo.length())))) { + group = it->first; + forwardSeq.setUnaligned(rawFSequence.substr(foligo.length())); + reverseSeq.setUnaligned(rawRSequence.substr(0,(rawRSequence.length()-roligo.length()))); + forwardQual.trimQScores(foligo.length(), -1); + reverseQual.trimQScores(-1, rawRSequence.length()-roligo.length()); + success = 0; + break; + } + } + + //if you found the barcode or if you don't want to allow for diffs + if ((bdiffs == 0) || (success == 0)) { return success; } + else { //try aligning and see if you can find it + + //look for forward + int maxLength = 0; + + Alignment* alignment; + if (ibarcodes.size() > 0) { + for(map::iterator it=ibarcodes.begin();it!=ibarcodes.end();it++){ + if(it->second.forward.length() > maxLength){ maxLength = it->second.forward.length(); } + } + alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); + }else{ alignment = NULL; } + + //can you find the barcode + int minDiff = 1e6; + int minCount = 1; + int minFGroup = -1; + int minFPos = 0; + + for(map::iterator it=ibarcodes.begin();it!=ibarcodes.end();it++){ + string oligo = it->second.forward; + + if(rawFSequence.length() < maxLength){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; + break; + } + + //use needleman to align first barcode.length()+numdiffs of sequence to each barcode + alignment->align(oligo, rawFSequence.substr(0,oligo.length()+bdiffs)); + oligo = alignment->getSeqAAln(); + string temp = alignment->getSeqBAln(); + + int alnLength = oligo.length(); + + for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } } + oligo = oligo.substr(0,alnLength); + temp = temp.substr(0,alnLength); + int numDiff = countDiffs(oligo, temp); + + if(numDiff < minDiff){ + minDiff = numDiff; + minCount = 1; + minFGroup = it->first; + minFPos = 0; + for(int i=0;i bdiffs) { success = minDiff; } //no good matches + else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes + else{ + //check for reverse match + if (alignment != NULL) { delete alignment; } + maxLength = 0; + + if (ibarcodes.size() > 0) { + for(map::iterator it=ibarcodes.begin();it!=ibarcodes.end();it++){ + if(it->second.reverse.length() > maxLength){ maxLength = it->second.reverse.length(); } + } + alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); + }else{ alignment = NULL; } + + //can you find the barcode + minDiff = 1e6; + minCount = 1; + int minRGroup = -1; + int minRPos = 0; + + for(map::iterator it=ibarcodes.begin();it!=ibarcodes.end();it++){ + string oligo = it->second.reverse; + + if(rawRSequence.length() < maxLength){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; + break; + } + + //use needleman to align first barcode.length()+numdiffs of sequence to each barcode + alignment->align(oligo, rawRSequence.substr((rawRSequence.length()-(oligo.length()+bdiffs)),oligo.length()+bdiffs)); + oligo = alignment->getSeqAAln(); + string temp = alignment->getSeqBAln(); + + int alnLength = oligo.length(); + for(int i=0;ifirst; + minRPos = 0; + for(int i=0;i bdiffs) { success = minDiff; } //no good matches + else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes + else{ + //we have an acceptable match for the forward and reverse, but do they match? + if (minFGroup == minRGroup) { + group = minFGroup; + forwardSeq.setUnaligned(rawFSequence.substr(minFPos)); + reverseSeq.setUnaligned(rawRSequence.substr(0,(rawRSequence.length()-minRPos))); + forwardQual.trimQScores(minFPos, -1); + reverseQual.trimQScores(-1, rawRSequence.length()-minRPos); + success = minDiff; + }else { success = bdiffs + 100; } + } + } + + if (alignment != NULL) { delete alignment; } + } + + return success; + + } + catch(exception& e) { + m->errorOut(e, "TrimOligos", "stripIBarcode"); + exit(1); + } } //*******************************************************************/ @@ -237,7 +428,7 @@ int TrimOligos::stripBarcode(Sequence& seq, int& group){ alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs)); oligo = alignment->getSeqAAln(); string temp = alignment->getSeqBAln(); - + int alnLength = oligo.length(); for(int i=oligo.length()-1;i>=0;i--){ @@ -245,7 +436,7 @@ int TrimOligos::stripBarcode(Sequence& seq, int& group){ } oligo = oligo.substr(0,alnLength); temp = temp.substr(0,alnLength); - + int numDiff = countDiffs(oligo, temp); if(numDiff < minDiff){ @@ -285,6 +476,239 @@ int TrimOligos::stripBarcode(Sequence& seq, int& group){ exit(1); } +} +/******************************************************************* +int TrimOligos::stripRBarcode(Sequence& seq, QualityScores& qual, int& group){ + try { + + string rawSequence = seq.getUnaligned(); + int success = bdiffs + 1; //guilty until proven innocent + + //can you find the barcode + for(map::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){ + string oligo = it->first; + if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out + break; + } + + if(compareDNASeq(oligo, rawSequence.substr((rawSequence.length()-oligo.length()),oligo.length()))){ + group = it->second; + seq.setUnaligned(rawSequence.substr(0,(rawSequence.length()-oligo.length()))); + + if(qual.getName() != ""){ + qual.trimQScores(-1, rawSequence.length()-oligo.length()); + } + + success = 0; + break; + } + } + + //if you found the barcode or if you don't want to allow for diffs + if ((bdiffs == 0) || (success == 0)) { return success; } + + else { //try aligning and see if you can find it + + int maxLength = 0; + + Alignment* alignment; + if (rbarcodes.size() > 0) { + map::iterator it; + + for(it=rbarcodes.begin();it!=rbarcodes.end();it++){ + if(it->first.length() > maxLength){ + maxLength = it->first.length(); + } + } + alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); + + }else{ alignment = NULL; } + + //can you find the barcode + int minDiff = 1e6; + int minCount = 1; + int minGroup = -1; + int minPos = 0; + + for(map::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){ + string oligo = it->first; + // int length = oligo.length(); + + if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; + break; + } + + //use needleman to align first barcode.length()+numdiffs of sequence to each barcode + alignment->align(oligo, rawSequence.substr((rawSequence.length()-(oligo.length()+bdiffs)),oligo.length()+bdiffs)); + oligo = alignment->getSeqAAln(); + string temp = alignment->getSeqBAln(); + + int alnLength = oligo.length(); + + for(int i=0;isecond; + minPos = 0; + for(int i=alnLength-1;i>=0;i--){ + if(temp[i] != '-'){ + minPos++; + } + } + } + else if(numDiff == minDiff){ + minCount++; + } + + } + + if(minDiff > bdiffs) { success = minDiff; } //no good matches + else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes + else{ //use the best match + group = minGroup; + seq.setUnaligned(rawSequence.substr(0, (rawSequence.length()-minPos))); + + if(qual.getName() != ""){ + qual.trimQScores(-1, (rawSequence.length()-minPos)); + } + success = minDiff; + } + + if (alignment != NULL) { delete alignment; } + + } + + return success; + + } + catch(exception& e) { + m->errorOut(e, "TrimOligos", "stripRBarcode"); + exit(1); + } + +} +/******************************************************************* +int TrimOligos::stripRBarcode(Sequence& seq, int& group){ + try { + + string rawSequence = seq.getUnaligned(); + int success = bdiffs + 1; //guilty until proven innocent + + //can you find the barcode + for(map::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){ + string oligo = it->first; + if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out + break; + } + + if(compareDNASeq(oligo, rawSequence.substr((rawSequence.length()-oligo.length()),oligo.length()))){ + group = it->second; + seq.setUnaligned(rawSequence.substr(0,(rawSequence.length()-oligo.length()))); + + success = 0; + break; + } + } + + //if you found the barcode or if you don't want to allow for diffs + if ((bdiffs == 0) || (success == 0)) { return success; } + + else { //try aligning and see if you can find it + + int maxLength = 0; + + Alignment* alignment; + if (rbarcodes.size() > 0) { + map::iterator it; + + for(it=rbarcodes.begin();it!=rbarcodes.end();it++){ + if(it->first.length() > maxLength){ + maxLength = it->first.length(); + } + } + alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1)); + + }else{ alignment = NULL; } + + //can you find the barcode + int minDiff = 1e6; + int minCount = 1; + int minGroup = -1; + int minPos = 0; + + for(map::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){ + string oligo = it->first; + // int length = oligo.length(); + + if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length + success = bdiffs + 10; + break; + } + + //use needleman to align first barcode.length()+numdiffs of sequence to each barcode + alignment->align(oligo, rawSequence.substr((rawSequence.length()-(oligo.length()+bdiffs)),oligo.length()+bdiffs)); + oligo = alignment->getSeqAAln(); + string temp = alignment->getSeqBAln(); + + int alnLength = oligo.length(); + + for(int i=0;isecond; + minPos = 0; + for(int i=alnLength-1;i>=0;i--){ + if(temp[i] != '-'){ + minPos++; + } + } + } + else if(numDiff == minDiff){ + minCount++; + } + + } + + if(minDiff > bdiffs) { success = minDiff; } //no good matches + else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes + else{ //use the best match + group = minGroup; + seq.setUnaligned(rawSequence.substr(0, (rawSequence.length()-minPos))); + + success = minDiff; + } + + if (alignment != NULL) { delete alignment; } + + } + + return success; + + } + catch(exception& e) { + m->errorOut(e, "TrimOligos", "stripRBarcode"); + exit(1); + } + } //********************************************************************/ int TrimOligos::stripForward(Sequence& seq, int& group){ @@ -576,13 +1000,13 @@ bool TrimOligos::stripReverse(Sequence& seq){ bool TrimOligos::stripLinker(Sequence& seq, QualityScores& qual){ try { string rawSequence = seq.getUnaligned(); - bool success = 0; //guilty until proven innocent + bool success = ldiffs + 1; //guilty until proven innocent for(int i=0;i