/**************************************************************************************************/
-NeedlemanOverlap::NeedlemanOverlap(float gO, float m, float mm, int r) :// note that we don't have a gap extend
-gap(gO), match(m), mismatch(mm), Alignment(r) { // the gap openning penalty is assessed for
- // every gapped position
- for(int i=1;i<nCols;i++){
- alignment[0][i].prevCell = 'l'; // initialize first row by pointing all poiters to the left
- alignment[0][i].cValue = 0; // and the score to zero
- }
+NeedlemanOverlap::NeedlemanOverlap(float gO, float f, float mm, int r) :// note that we don't have a gap extend
+gap(gO), match(f), mismatch(mm), Alignment(r) { // the gap openning penalty is assessed for
+ try { // every gapped position
+ for(int i=1;i<nCols;i++){
+ alignment[0][i].prevCell = 'l'; // initialize first row by pointing all poiters to the left
+ alignment[0][i].cValue = 0; // and the score to zero
+ }
+
+ for(int i=1;i<nRows;i++){
+ alignment[i][0].prevCell = 'u'; // initialize first column by pointing all poiters upwards
+ alignment[i][0].cValue = 0; // and the score to zero
+ }
- for(int i=1;i<nRows;i++){
- alignment[i][0].prevCell = 'u'; // initialize first column by pointing all poiters upwards
- alignment[i][0].cValue = 0; // and the score to zero
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NeedlemanOverlap", "NeedlemanOverlap");
+ exit(1);
}
}
-
/**************************************************************************************************/
NeedlemanOverlap::~NeedlemanOverlap(){ /* do nothing */ }
/**************************************************************************************************/
void NeedlemanOverlap::align(string A, string B){
+ try {
- seqA = ' ' + A; lA = seqA.length(); // algorithm requires a dummy space at the beginning of each string
- seqB = ' ' + B; lB = seqB.length(); // algorithm requires a dummy space at the beginning of each string
+ seqA = ' ' + A; lA = seqA.length(); // algorithm requires a dummy space at the beginning of each string
+ seqB = ' ' + B; lB = seqB.length(); // algorithm requires a dummy space at the beginning of each string
+
+ if (lA > nRows) { m->mothurOut("One of your candidate sequences is longer than you longest template sequence. Your longest template sequence is " + toString(nRows) + ". Your candidate is " + toString(lA) + "."); m->mothurOutEndLine(); }
+
+ for(int i=1;i<lB;i++){ // This code was largely translated from Perl code provided in Ex 3.1
+
+ for(int j=1;j<lA;j++){ // of the O'Reilly BLAST book. I found that the example output had a
- for(int i=1;i<lB;i++){ // This code was largely translated from Perl code provided in Ex 3.1
- for(int j=1;j<lA;j++){ // of the O'Reilly BLAST book. I found that the example output had a
- // number of errors
- float diagonal;
- if(seqB[i] == seqA[j]) { diagonal = alignment[i-1][j-1].cValue + match; }
- else { diagonal = alignment[i-1][j-1].cValue + mismatch; }
-
- float up = alignment[i-1][j].cValue + gap;
- float left = alignment[i][j-1].cValue + gap;
+ // number of errors
+ float diagonal;
+ if(seqB[i] == seqA[j]) { diagonal = alignment[i-1][j-1].cValue + match; }
+ else { diagonal = alignment[i-1][j-1].cValue + mismatch; }
- if(diagonal >= up){
- if(diagonal >= left){
- alignment[i][j].cValue = diagonal;
- alignment[i][j].prevCell = 'd';
+ float up = alignment[i-1][j].cValue + gap;
+ float left = alignment[i][j-1].cValue + gap;
+
+ if(diagonal >= up){
+ if(diagonal >= left){
+ alignment[i][j].cValue = diagonal;
+ alignment[i][j].prevCell = 'd';
+ }
+ else{
+ alignment[i][j].cValue = left;
+ alignment[i][j].prevCell = 'l';
+ }
}
else{
- alignment[i][j].cValue = left;
- alignment[i][j].prevCell = 'l';
+ if(up >= left){
+ alignment[i][j].cValue = up;
+ alignment[i][j].prevCell = 'u';
+ }
+ else{
+ alignment[i][j].cValue = left;
+ alignment[i][j].prevCell = 'l';
+ }
}
}
- else{
- if(up >= left){
- alignment[i][j].cValue = up;
- alignment[i][j].prevCell = 'u';
+ }
+
+ Overlap over;
+ over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences
+ traceBack(); // Traceback the alignment to populate seqAaln and seqBaln
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NeedlemanOverlap", "align");
+ exit(1);
+ }
+
+}
+/**************************************************************************************************/
+
+void NeedlemanOverlap::alignPrimer(string A, string B){
+ try {
+
+ seqA = ' ' + A; lA = seqA.length(); // algorithm requires a dummy space at the beginning of each string
+ seqB = ' ' + B; lB = seqB.length(); // algorithm requires a dummy space at the beginning of each string
+
+ if (lA > nRows) { m->mothurOut("One of your candidate sequences is longer than you longest template sequence. Your longest template sequence is " + toString(nRows) + ". Your candidate is " + toString(lA) + "."); m->mothurOutEndLine(); }
+
+ for(int i=1;i<lB;i++){ // This code was largely translated from Perl code provided in Ex 3.1
+
+ for(int j=1;j<lA;j++){ // of the O'Reilly BLAST book. I found that the example output had a
+
+ // number of errors
+ float diagonal;
+ if(isEquivalent(seqB[i],seqA[j])) { diagonal = alignment[i-1][j-1].cValue + match; }
+ else { diagonal = alignment[i-1][j-1].cValue + mismatch; }
+
+ float up = alignment[i-1][j].cValue + gap;
+ float left = alignment[i][j-1].cValue + gap;
+
+ if(diagonal >= up){
+ if(diagonal >= left){
+ alignment[i][j].cValue = diagonal;
+ alignment[i][j].prevCell = 'd';
+ }
+ else{
+ alignment[i][j].cValue = left;
+ alignment[i][j].prevCell = 'l';
+ }
}
else{
- alignment[i][j].cValue = left;
- alignment[i][j].prevCell = 'l';
+ if(up >= left){
+ alignment[i][j].cValue = up;
+ alignment[i][j].prevCell = 'u';
+ }
+ else{
+ alignment[i][j].cValue = left;
+ alignment[i][j].prevCell = 'l';
+ }
}
}
}
+
+ Overlap over;
+ over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences
+ traceBack(); // Traceback the alignment to populate seqAaln and seqBaln
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "NeedlemanOverlap", "alignPrimer");
+ exit(1);
}
- Overlap over;
- over.setOverlap(alignment, lA, lB, 0); // Fix gaps at the beginning and end of the sequences
- traceBack(); // Traceback the alignment to populate seqAaln and seqBaln
+
}
+//********************************************************************/
+bool NeedlemanOverlap::isEquivalent(char oligo, char seq){
+ try {
+
+ bool same = true;
+
+ oligo = toupper(oligo);
+ seq = toupper(seq);
+
+ if(oligo != seq){
+ if(oligo == 'A' && (seq != 'A' && seq != 'M' && seq != 'R' && seq != 'W' && seq != 'D' && seq != 'H' && seq != 'V')) { same = false; }
+ else if(oligo == 'C' && (seq != 'C' && seq != 'Y' && seq != 'M' && seq != 'S' && seq != 'B' && seq != 'H' && seq != 'V')) { same = false; }
+ else if(oligo == 'G' && (seq != 'G' && seq != 'R' && seq != 'K' && seq != 'S' && seq != 'B' && seq != 'D' && seq != 'V')) { same = false; }
+ else if(oligo == 'T' && (seq != 'T' && seq != 'Y' && seq != 'K' && seq != 'W' && seq != 'B' && seq != 'D' && seq != 'H')) { same = false; }
+ else if((oligo == '.' || oligo == '-')) { same = false; }
+ else if((oligo == 'N' || oligo == 'I') && (seq == 'N')) { same = false; }
+ else if(oligo == 'R' && (seq != 'A' && seq != 'G')) { same = false; }
+ else if(oligo == 'Y' && (seq != 'C' && seq != 'T')) { same = false; }
+ else if(oligo == 'M' && (seq != 'C' && seq != 'A')) { same = false; }
+ else if(oligo == 'K' && (seq != 'T' && seq != 'G')) { same = false; }
+ else if(oligo == 'W' && (seq != 'T' && seq != 'A')) { same = false; }
+ else if(oligo == 'S' && (seq != 'C' && seq != 'G')) { same = false; }
+ else if(oligo == 'B' && (seq != 'C' && seq != 'T' && seq != 'G')) { same = false; }
+ else if(oligo == 'D' && (seq != 'A' && seq != 'T' && seq != 'G')) { same = false; }
+ else if(oligo == 'H' && (seq != 'A' && seq != 'T' && seq != 'C')) { same = false; }
+ else if(oligo == 'V' && (seq != 'A' && seq != 'C' && seq != 'G')) { same = false; }
+ }
+
+
+ return same;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "TrimOligos", "countDiffs");
+ exit(1);
+ }
+}
/**************************************************************************************************/
projects / mothur.git / blobdiff