trimoligos.cpp

   1 /*
   2  *  trimoligos.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 9/1/11.
   6  *  Copyright 2011 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "trimoligos.h"
  11 #include "alignment.hpp"
  12 #include "needlemanoverlap.hpp"
  13
  14
  15 /********************************************************************/
  16 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
  17 TrimOligos::TrimOligos(int p, int b, int l, int s, map<string, int> pr, map<string, int> br, map<string, int> rbr, vector<string> r, vector<string> lk, vector<string> sp){
  18         try {
  19                 m = MothurOut::getInstance();
  20
  21                 pdiffs = p;
  22                 bdiffs = b;
  23         ldiffs = l;
  24         sdiffs = s;
  25
  26                 barcodes = br;
  27         rbarcodes = rbr;
  28                 primers = pr;
  29                 revPrimer = r;
  30         linker = lk;
  31         spacer = sp;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "TrimOligos", "TrimOligos");
  35                 exit(1);
  36         }
  37 }
  38 /********************************************************************/
  39 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
  40 TrimOligos::TrimOligos(int p, int b, int l, int s, map<string, int> pr, map<string, int> br, vector<string> r, vector<string> lk, vector<string> sp){
  41         try {
  42                 m = MothurOut::getInstance();
  43
  44                 pdiffs = p;
  45                 bdiffs = b;
  46         ldiffs = l;
  47         sdiffs = s;
  48
  49                 barcodes = br;
  50                 primers = pr;
  51                 revPrimer = r;
  52         linker = lk;
  53         spacer = sp;
  54         }
  55         catch(exception& e) {
  56                 m->errorOut(e, "TrimOligos", "TrimOligos");
  57                 exit(1);
  58         }
  59 }
  60 /********************************************************************/
  61 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
  62 TrimOligos::TrimOligos(int p, int b, map<string, int> pr, map<string, int> br, vector<string> r){
  63         try {
  64                 m = MothurOut::getInstance();
  65
  66                 pdiffs = p;
  67                 bdiffs = b;
  68
  69                 barcodes = br;
  70                 primers = pr;
  71                 revPrimer = r;
  72         }
  73         catch(exception& e) {
  74                 m->errorOut(e, "TrimOligos", "TrimOligos");
  75                 exit(1);
  76         }
  77 }
  78 /********************************************************************/
  79 TrimOligos::~TrimOligos() {}
  80 //*******************************************************************/
  81 int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){
  82         try {
  83
  84                 string rawSequence = seq.getUnaligned();
  85                 int success = bdiffs + 1;       //guilty until proven innocent
  86
  87                 //can you find the barcode
  88                 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
  89                         string oligo = it->first;
  90                         if(rawSequence.length() < oligo.length()){      //let's just assume that the barcodes are the same length
  91                                 success = bdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
  92                                 break;
  93                         }
  94
  95                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
  96                                 group = it->second;
  97                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
  98
  99                                 if(qual.getName() != ""){
 100                                         qual.trimQScores(oligo.length(), -1);
 101                                 }
 102
 103                                 success = 0;
 104                                 break;
 105                         }
 106                 }
 107
 108                 //if you found the barcode or if you don't want to allow for diffs
 109                 if ((bdiffs == 0) || (success == 0)) { return success;  }
 110
 111                 else { //try aligning and see if you can find it
 112
 113                         int maxLength = 0;
 114
 115                         Alignment* alignment;
 116                         if (barcodes.size() > 0) {
 117                                 map<string,int>::iterator it;
 118
 119                                 for(it=barcodes.begin();it!=barcodes.end();it++){
 120                                         if(it->first.length() > maxLength){
 121                                                 maxLength = it->first.length();
 122                                         }
 123                                 }
 124                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
 125
 126                         }else{ alignment = NULL; }
 127
 128                         //can you find the barcode
 129                         int minDiff = 1e6;
 130                         int minCount = 1;
 131                         int minGroup = -1;
 132                         int minPos = 0;
 133
 134                         for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 135                                 string oligo = it->first;
 136                                 //                              int length = oligo.length();
 137
 138                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 139                                         success = bdiffs + 10;
 140                                         break;
 141                                 }
 142
 143                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 144                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
 145                                 oligo = alignment->getSeqAAln();
 146                                 string temp = alignment->getSeqBAln();
 147
 148                                 int alnLength = oligo.length();
 149
 150                                 for(int i=oligo.length()-1;i>=0;i--){
 151                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 152                                 }
 153                                 oligo = oligo.substr(0,alnLength);
 154                                 temp = temp.substr(0,alnLength);
 155                 int numDiff = countDiffs(oligo, temp);
 156
 157                                 if(numDiff < minDiff){
 158                                         minDiff = numDiff;
 159                                         minCount = 1;
 160                                         minGroup = it->second;
 161                                         minPos = 0;
 162                                         for(int i=0;i<alnLength;i++){
 163                                                 if(temp[i] != '-'){
 164                                                         minPos++;
 165                                                 }
 166                                         }
 167                                 }
 168                                 else if(numDiff == minDiff){
 169                                         minCount++;
 170                                 }
 171
 172                         }
 173
 174                         if(minDiff > bdiffs)    {       success = minDiff;              }       //no good matches
 175                         else if(minCount > 1)   {       success = bdiffs + 100; }       //can't tell the difference between multiple barcodes
 176                         else{                                                                                                   //use the best match
 177                                 group = minGroup;
 178                                 seq.setUnaligned(rawSequence.substr(minPos));
 179
 180                                 if(qual.getName() != ""){
 181                                         qual.trimQScores(minPos, -1);
 182                                 }
 183                                 success = minDiff;
 184                         }
 185
 186                         if (alignment != NULL) {  delete alignment;  }
 187
 188                 }
 189
 190                 return success;
 191
 192         }
 193         catch(exception& e) {
 194                 m->errorOut(e, "TrimOligos", "stripBarcode");
 195                 exit(1);
 196         }
 197
 198 }
 199 //*******************************************************************/
 200 int TrimOligos::stripBarcode(Sequence& seq, int& group){
 201         try {
 202
 203                 string rawSequence = seq.getUnaligned();
 204                 int success = bdiffs + 1;       //guilty until proven innocent
 205
 206                 //can you find the barcode
 207                 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 208                         string oligo = it->first;
 209                         if(rawSequence.length() < oligo.length()){      //let's just assume that the barcodes are the same length
 210                                 success = bdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 211                                 break;
 212                         }
 213
 214                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 215                                 group = it->second;
 216                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 217
 218                                 success = 0;
 219                                 break;
 220                         }
 221                 }
 222
 223                 //if you found the barcode or if you don't want to allow for diffs
 224                 if ((bdiffs == 0) || (success == 0)) { return success;  }
 225
 226                 else { //try aligning and see if you can find it
 227
 228                         int maxLength = 0;
 229
 230                         Alignment* alignment;
 231                         if (barcodes.size() > 0) {
 232                                 map<string,int>::iterator it=barcodes.begin();
 233
 234                                 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 235                                         if(it->first.length() > maxLength){
 236                                                 maxLength = it->first.length();
 237                                         }
 238                                 }
 239                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
 240
 241                         }else{ alignment = NULL; }
 242
 243                         //can you find the barcode
 244                         int minDiff = 1e6;
 245                         int minCount = 1;
 246                         int minGroup = -1;
 247                         int minPos = 0;
 248
 249                         for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 250                                 string oligo = it->first;
 251                                 //                              int length = oligo.length();
 252
 253                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 254                                         success = bdiffs + 10;
 255                                         break;
 256                                 }
 257
 258                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 259                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
 260                                 oligo = alignment->getSeqAAln();
 261                                 string temp = alignment->getSeqBAln();
 262
 263                                 int alnLength = oligo.length();
 264
 265                                 for(int i=oligo.length()-1;i>=0;i--){
 266                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 267                                 }
 268                                 oligo = oligo.substr(0,alnLength);
 269                                 temp = temp.substr(0,alnLength);
 270
 271                                 int numDiff = countDiffs(oligo, temp);
 272
 273                                 if(numDiff < minDiff){
 274                                         minDiff = numDiff;
 275                                         minCount = 1;
 276                                         minGroup = it->second;
 277                                         minPos = 0;
 278                                         for(int i=0;i<alnLength;i++){
 279                                                 if(temp[i] != '-'){
 280                                                         minPos++;
 281                                                 }
 282                                         }
 283                                 }
 284                                 else if(numDiff == minDiff){
 285                                         minCount++;
 286                                 }
 287
 288                         }
 289
 290                         if(minDiff > bdiffs)    {       success = minDiff;              }       //no good matches
 291                         else if(minCount > 1)   {       success = bdiffs + 100; }       //can't tell the difference between multiple barcodes
 292                         else{                                                                                                   //use the best match
 293                                 group = minGroup;
 294                                 seq.setUnaligned(rawSequence.substr(minPos));
 295                                 success = minDiff;
 296                         }
 297
 298                         if (alignment != NULL) {  delete alignment;  }
 299
 300                 }
 301
 302                 return success;
 303
 304         }
 305         catch(exception& e) {
 306                 m->errorOut(e, "TrimOligos", "stripBarcode");
 307                 exit(1);
 308         }
 309
 310 }
 311 //*******************************************************************/
 312 int TrimOligos::stripRBarcode(Sequence& seq, QualityScores& qual, int& group){
 313         try {
 314
 315                 string rawSequence = seq.getUnaligned();
 316                 int success = bdiffs + 1;       //guilty until proven innocent
 317
 318                 //can you find the barcode
 319                 for(map<string,int>::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){
 320                         string oligo = it->first;
 321                         if(rawSequence.length() < oligo.length()){      //let's just assume that the barcodes are the same length
 322                                 success = bdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 323                                 break;
 324                         }
 325
 326                         if(compareDNASeq(oligo, rawSequence.substr((rawSequence.length()-oligo.length()),oligo.length()))){
 327                                 group = it->second;
 328                                 seq.setUnaligned(rawSequence.substr(0,(rawSequence.length()-oligo.length())));
 329
 330                                 if(qual.getName() != ""){
 331                                         qual.trimQScores(-1, rawSequence.length()-oligo.length());
 332                                 }
 333
 334                                 success = 0;
 335                                 break;
 336                         }
 337                 }
 338
 339                 //if you found the barcode or if you don't want to allow for diffs
 340                 if ((bdiffs == 0) || (success == 0)) { return success;  }
 341
 342                 else { //try aligning and see if you can find it
 343
 344                         int maxLength = 0;
 345
 346                         Alignment* alignment;
 347                         if (rbarcodes.size() > 0) {
 348                                 map<string,int>::iterator it;
 349
 350                                 for(it=rbarcodes.begin();it!=rbarcodes.end();it++){
 351                                         if(it->first.length() > maxLength){
 352                                                 maxLength = it->first.length();
 353                                         }
 354                                 }
 355                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
 356
 357                         }else{ alignment = NULL; }
 358
 359                         //can you find the barcode
 360                         int minDiff = 1e6;
 361                         int minCount = 1;
 362                         int minGroup = -1;
 363                         int minPos = 0;
 364
 365                         for(map<string,int>::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){
 366                                 string oligo = it->first;
 367                                 //                              int length = oligo.length();
 368
 369                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 370                                         success = bdiffs + 10;
 371                                         break;
 372                                 }
 373
 374                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 375                                 alignment->align(oligo, rawSequence.substr((rawSequence.length()-(oligo.length()+bdiffs)),oligo.length()+bdiffs));
 376                                 oligo = alignment->getSeqAAln();
 377                                 string temp = alignment->getSeqBAln();
 378
 379                                 int alnLength = oligo.length();
 380
 381                                 for(int i=0;i<alnLength;i++){
 382                                         if(oligo[i] != '-'){    alnLength = i;  break;  }
 383                                 }
 384                                 oligo = oligo.substr(alnLength);
 385                                 temp = temp.substr(alnLength);
 386
 387                                 int numDiff = countDiffs(oligo, temp);
 388
 389                                 if(numDiff < minDiff){
 390                                         minDiff = numDiff;
 391                                         minCount = 1;
 392                                         minGroup = it->second;
 393                                         minPos = 0;
 394                                         for(int i=alnLength-1;i>=0;i--){
 395                                                 if(temp[i] != '-'){
 396                                                         minPos++;
 397                                                 }
 398                                         }
 399                                 }
 400                                 else if(numDiff == minDiff){
 401                                         minCount++;
 402                                 }
 403
 404                         }
 405
 406                         if(minDiff > bdiffs)    {       success = minDiff;              }       //no good matches
 407                         else if(minCount > 1)   {       success = bdiffs + 100; }       //can't tell the difference between multiple barcodes
 408                         else{                                                                                                   //use the best match
 409                                 group = minGroup;
 410                                 seq.setUnaligned(rawSequence.substr(0, (rawSequence.length()-minPos)));
 411
 412                                 if(qual.getName() != ""){
 413                                         qual.trimQScores(-1, (rawSequence.length()-minPos));
 414                                 }
 415                                 success = minDiff;
 416                         }
 417
 418                         if (alignment != NULL) {  delete alignment;  }
 419
 420                 }
 421
 422                 return success;
 423
 424         }
 425         catch(exception& e) {
 426                 m->errorOut(e, "TrimOligos", "stripRBarcode");
 427                 exit(1);
 428         }
 429
 430 }
 431 //*******************************************************************/
 432 int TrimOligos::stripRBarcode(Sequence& seq, int& group){
 433         try {
 434
 435                 string rawSequence = seq.getUnaligned();
 436                 int success = bdiffs + 1;       //guilty until proven innocent
 437
 438                 //can you find the barcode
 439                 for(map<string,int>::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){
 440                         string oligo = it->first;
 441                         if(rawSequence.length() < oligo.length()){      //let's just assume that the barcodes are the same length
 442                                 success = bdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 443                                 break;
 444                         }
 445
 446                         if(compareDNASeq(oligo, rawSequence.substr((rawSequence.length()-oligo.length()),oligo.length()))){
 447                                 group = it->second;
 448                                 seq.setUnaligned(rawSequence.substr(0,(rawSequence.length()-oligo.length())));
 449
 450                                 success = 0;
 451                                 break;
 452                         }
 453                 }
 454
 455                 //if you found the barcode or if you don't want to allow for diffs
 456                 if ((bdiffs == 0) || (success == 0)) { return success;  }
 457
 458                 else { //try aligning and see if you can find it
 459
 460                         int maxLength = 0;
 461
 462                         Alignment* alignment;
 463                         if (rbarcodes.size() > 0) {
 464                                 map<string,int>::iterator it;
 465
 466                                 for(it=rbarcodes.begin();it!=rbarcodes.end();it++){
 467                                         if(it->first.length() > maxLength){
 468                                                 maxLength = it->first.length();
 469                                         }
 470                                 }
 471                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
 472
 473                         }else{ alignment = NULL; }
 474
 475                         //can you find the barcode
 476                         int minDiff = 1e6;
 477                         int minCount = 1;
 478                         int minGroup = -1;
 479                         int minPos = 0;
 480
 481                         for(map<string,int>::iterator it=rbarcodes.begin();it!=rbarcodes.end();it++){
 482                                 string oligo = it->first;
 483                                 //                              int length = oligo.length();
 484
 485                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 486                                         success = bdiffs + 10;
 487                                         break;
 488                                 }
 489
 490                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 491                                 alignment->align(oligo, rawSequence.substr((rawSequence.length()-(oligo.length()+bdiffs)),oligo.length()+bdiffs));
 492                                 oligo = alignment->getSeqAAln();
 493                                 string temp = alignment->getSeqBAln();
 494
 495                                 int alnLength = oligo.length();
 496
 497                                 for(int i=0;i<alnLength;i++){
 498                                         if(oligo[i] != '-'){    alnLength = i;  break;  }
 499                                 }
 500                                 oligo = oligo.substr(alnLength);
 501                                 temp = temp.substr(alnLength);
 502
 503                                 int numDiff = countDiffs(oligo, temp);
 504
 505                                 if(numDiff < minDiff){
 506                                         minDiff = numDiff;
 507                                         minCount = 1;
 508                                         minGroup = it->second;
 509                                         minPos = 0;
 510                                         for(int i=alnLength-1;i>=0;i--){
 511                                                 if(temp[i] != '-'){
 512                                                         minPos++;
 513                                                 }
 514                                         }
 515                                 }
 516                                 else if(numDiff == minDiff){
 517                                         minCount++;
 518                                 }
 519
 520                         }
 521
 522                         if(minDiff > bdiffs)    {       success = minDiff;              }       //no good matches
 523                         else if(minCount > 1)   {       success = bdiffs + 100; }       //can't tell the difference between multiple barcodes
 524                         else{                                                                                                   //use the best match
 525                                 group = minGroup;
 526                                 seq.setUnaligned(rawSequence.substr(0, (rawSequence.length()-minPos)));
 527
 528                                 success = minDiff;
 529                         }
 530
 531                         if (alignment != NULL) {  delete alignment;  }
 532
 533                 }
 534
 535                 return success;
 536
 537         }
 538         catch(exception& e) {
 539                 m->errorOut(e, "TrimOligos", "stripRBarcode");
 540                 exit(1);
 541         }
 542
 543 }
 544 //********************************************************************/
 545 int TrimOligos::stripForward(Sequence& seq, int& group){
 546         try {
 547
 548                 string rawSequence = seq.getUnaligned();
 549                 int success = pdiffs + 1;       //guilty until proven innocent
 550
 551                 //can you find the primer
 552                 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 553                         string oligo = it->first;
 554                         if(rawSequence.length() < oligo.length()){      //let's just assume that the primers are the same length
 555                                 success = pdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 556                                 break;
 557                         }
 558
 559                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 560                                 group = it->second;
 561                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 562                                 success = 0;
 563                                 break;
 564                         }
 565                 }
 566
 567                 //if you found the barcode or if you don't want to allow for diffs
 568                 if ((pdiffs == 0) || (success == 0)) {  return success;  }
 569
 570                 else { //try aligning and see if you can find it
 571
 572                         int maxLength = 0;
 573
 574                         Alignment* alignment;
 575                         if (primers.size() > 0) {
 576                                 map<string,int>::iterator it;
 577
 578                                 for(it=primers.begin();it!=primers.end();it++){
 579                                         if(it->first.length() > maxLength){
 580                                                 maxLength = it->first.length();
 581                                         }
 582                                 }
 583                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
 584
 585                         }else{ alignment = NULL; }
 586
 587                         //can you find the barcode
 588                         int minDiff = 1e6;
 589                         int minCount = 1;
 590                         int minGroup = -1;
 591                         int minPos = 0;
 592
 593                         for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 594                                 string oligo = it->first;
 595                                 //                              int length = oligo.length();
 596
 597                                 if(rawSequence.length() < maxLength){
 598                                         success = pdiffs + 100;
 599                                         break;
 600                                 }
 601
 602                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 603                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
 604                                 oligo = alignment->getSeqAAln();
 605                                 string temp = alignment->getSeqBAln();
 606
 607                                 int alnLength = oligo.length();
 608
 609                                 for(int i=oligo.length()-1;i>=0;i--){
 610                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 611                                 }
 612                                 oligo = oligo.substr(0,alnLength);
 613                                 temp = temp.substr(0,alnLength);
 614
 615                                 int numDiff = countDiffs(oligo, temp);
 616
 617                                 if(numDiff < minDiff){
 618                                         minDiff = numDiff;
 619                                         minCount = 1;
 620                                         minGroup = it->second;
 621                                         minPos = 0;
 622                                         for(int i=0;i<alnLength;i++){
 623                                                 if(temp[i] != '-'){
 624                                                         minPos++;
 625                                                 }
 626                                         }
 627                                 }
 628                                 else if(numDiff == minDiff){
 629                                         minCount++;
 630                                 }
 631
 632                         }
 633
 634                         if(minDiff > pdiffs)    {       success = minDiff;              }       //no good matches
 635                         else if(minCount > 1)   {       success = pdiffs + 10;  }       //can't tell the difference between multiple primers
 636                         else{                                                                                                   //use the best match
 637                                 group = minGroup;
 638                                 seq.setUnaligned(rawSequence.substr(minPos));
 639                                 success = minDiff;
 640                         }
 641
 642                         if (alignment != NULL) {  delete alignment;  }
 643
 644                 }
 645
 646                 return success;
 647
 648         }
 649         catch(exception& e) {
 650                 m->errorOut(e, "TrimOligos", "stripForward");
 651                 exit(1);
 652         }
 653 }
 654 //*******************************************************************/
 655 int TrimOligos::stripForward(Sequence& seq, QualityScores& qual, int& group, bool keepForward){
 656         try {
 657                 string rawSequence = seq.getUnaligned();
 658                 int success = pdiffs + 1;       //guilty until proven innocent
 659
 660                 //can you find the primer
 661                 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 662                         string oligo = it->first;
 663                         if(rawSequence.length() < oligo.length()){      //let's just assume that the primers are the same length
 664                                 success = pdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 665                                 break;
 666                         }
 667
 668                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 669                                 group = it->second;
 670                                 if (!keepForward) { seq.setUnaligned(rawSequence.substr(oligo.length())); }
 671                                 if(qual.getName() != ""){
 672                                         if (!keepForward) {  qual.trimQScores(oligo.length(), -1); }
 673                                 }
 674                                 success = 0;
 675                                 break;
 676                         }
 677                 }
 678
 679                 //if you found the barcode or if you don't want to allow for diffs
 680                 if ((pdiffs == 0) || (success == 0)) { return success;  }
 681
 682                 else { //try aligning and see if you can find it
 683
 684                         int maxLength = 0;
 685
 686                         Alignment* alignment;
 687                         if (primers.size() > 0) {
 688                                 map<string,int>::iterator it;
 689
 690                                 for(it=primers.begin();it!=primers.end();it++){
 691                                         if(it->first.length() > maxLength){
 692                                                 maxLength = it->first.length();
 693                                         }
 694                                 }
 695                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
 696
 697                         }else{ alignment = NULL; }
 698
 699                         //can you find the barcode
 700                         int minDiff = 1e6;
 701                         int minCount = 1;
 702                         int minGroup = -1;
 703                         int minPos = 0;
 704
 705                         for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 706                                 string oligo = it->first;
 707                                 //                              int length = oligo.length();
 708
 709                                 if(rawSequence.length() < maxLength){
 710                                         success = pdiffs + 100;
 711                                         break;
 712                                 }
 713
 714                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 715                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
 716                                 oligo = alignment->getSeqAAln();
 717                                 string temp = alignment->getSeqBAln();
 718
 719                                 int alnLength = oligo.length();
 720
 721                                 for(int i=oligo.length()-1;i>=0;i--){
 722                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 723                                 }
 724                                 oligo = oligo.substr(0,alnLength);
 725                                 temp = temp.substr(0,alnLength);
 726
 727                                 int numDiff = countDiffs(oligo, temp);
 728
 729                                 if(numDiff < minDiff){
 730                                         minDiff = numDiff;
 731                                         minCount = 1;
 732                                         minGroup = it->second;
 733                                         minPos = 0;
 734                                         for(int i=0;i<alnLength;i++){
 735                                                 if(temp[i] != '-'){
 736                                                         minPos++;
 737                                                 }
 738                                         }
 739                                 }
 740                                 else if(numDiff == minDiff){
 741                                         minCount++;
 742                                 }
 743
 744                         }
 745
 746                         if(minDiff > pdiffs)    {       success = minDiff;              }       //no good matches
 747                         else if(minCount > 1)   {       success = pdiffs + 10;  }       //can't tell the difference between multiple primers
 748                         else{                                                                                                   //use the best match
 749                                 group = minGroup;
 750                                 if (!keepForward) { seq.setUnaligned(rawSequence.substr(minPos)); }
 751                                 if(qual.getName() != ""){
 752                                         if (!keepForward) { qual.trimQScores(minPos, -1); }
 753                                 }
 754                                 success = minDiff;
 755                         }
 756
 757                         if (alignment != NULL) {  delete alignment;  }
 758
 759                 }
 760
 761                 return success;
 762
 763         }
 764         catch(exception& e) {
 765                 m->errorOut(e, "TrimOligos", "stripForward");
 766                 exit(1);
 767         }
 768 }
 769 //******************************************************************/
 770 bool TrimOligos::stripReverse(Sequence& seq, QualityScores& qual){
 771         try {
 772                 string rawSequence = seq.getUnaligned();
 773                 bool success = 0;       //guilty until proven innocent
 774
 775                 for(int i=0;i<revPrimer.size();i++){
 776                         string oligo = revPrimer[i];
 777
 778                         if(rawSequence.length() < oligo.length()){
 779                                 success = 0;
 780                                 break;
 781                         }
 782
 783                         if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
 784                                 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
 785                                 if(qual.getName() != ""){
 786                                         qual.trimQScores(-1, rawSequence.length()-oligo.length());
 787                                 }
 788                                 success = 1;
 789                                 break;
 790                         }
 791                 }
 792                 return success;
 793
 794         }
 795         catch(exception& e) {
 796                 m->errorOut(e, "TrimOligos", "stripReverse");
 797                 exit(1);
 798         }
 799 }
 800 //******************************************************************/
 801 bool TrimOligos::stripReverse(Sequence& seq){
 802         try {
 803
 804                 string rawSequence = seq.getUnaligned();
 805                 bool success = 0;       //guilty until proven innocent
 806
 807                 for(int i=0;i<revPrimer.size();i++){
 808                         string oligo = revPrimer[i];
 809
 810                         if(rawSequence.length() < oligo.length()){
 811                                 success = 0;
 812                                 break;
 813                         }
 814
 815                         if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
 816                                 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
 817                                 success = 1;
 818                                 break;
 819                         }
 820                 }
 821
 822                 return success;
 823
 824         }
 825         catch(exception& e) {
 826                 m->errorOut(e, "TrimOligos", "stripReverse");
 827                 exit(1);
 828         }
 829 }
 830 //******************************************************************/
 831 bool TrimOligos::stripLinker(Sequence& seq, QualityScores& qual){
 832         try {
 833                 string rawSequence = seq.getUnaligned();
 834                 bool success = ldiffs + 1;      //guilty until proven innocent
 835
 836                 for(int i=0;i<linker.size();i++){
 837                         string oligo = linker[i];
 838
 839                         if(rawSequence.length() < oligo.length()){
 840                                 success = ldiffs + 10;
 841                                 break;
 842                         }
 843
 844                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 845                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 846                                 if(qual.getName() != ""){
 847                                         qual.trimQScores(oligo.length(), -1);
 848                                 }
 849                                 success = 0;
 850                                 break;
 851                         }
 852                 }
 853
 854         //if you found the linker or if you don't want to allow for diffs
 855                 if ((ldiffs == 0) || (success == 0)) { return success;  }
 856
 857                 else { //try aligning and see if you can find it
 858
 859                         int maxLength = 0;
 860
 861                         Alignment* alignment;
 862                         if (linker.size() > 0) {
 863                                 for(int i = 0; i < linker.size(); i++){
 864                                         if(linker[i].length() > maxLength){
 865                                                 maxLength = linker[i].length();
 866                                         }
 867                                 }
 868                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+ldiffs+1));
 869
 870                         }else{ alignment = NULL; }
 871
 872                         //can you find the barcode
 873                         int minDiff = 1e6;
 874                         int minCount = 1;
 875                         int minPos = 0;
 876
 877                         for(int i = 0; i < linker.size(); i++){
 878                                 string oligo = linker[i];
 879                                 //                              int length = oligo.length();
 880
 881                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 882                                         success = ldiffs + 10;
 883                                         break;
 884                                 }
 885
 886                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 887                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
 888                                 oligo = alignment->getSeqAAln();
 889                                 string temp = alignment->getSeqBAln();
 890
 891                                 int alnLength = oligo.length();
 892
 893                                 for(int i=oligo.length()-1;i>=0;i--){
 894                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 895                                 }
 896                                 oligo = oligo.substr(0,alnLength);
 897                                 temp = temp.substr(0,alnLength);
 898
 899                                 int numDiff = countDiffs(oligo, temp);
 900
 901                                 if(numDiff < minDiff){
 902                                         minDiff = numDiff;
 903                                         minCount = 1;
 904                                         minPos = 0;
 905                                         for(int i=0;i<alnLength;i++){
 906                                                 if(temp[i] != '-'){
 907                                                         minPos++;
 908                                                 }
 909                                         }
 910                                 }
 911                                 else if(numDiff == minDiff){
 912                                         minCount++;
 913                                 }
 914
 915                         }
 916
 917                         if(minDiff > ldiffs)    {       success = minDiff;              }       //no good matches
 918                         else if(minCount > 1)   {       success = ldiffs + 100; }       //can't tell the difference between multiple barcodes
 919                         else{                                                                                                   //use the best match
 920                                 seq.setUnaligned(rawSequence.substr(minPos));
 921
 922                                 if(qual.getName() != ""){
 923                                         qual.trimQScores(minPos, -1);
 924                                 }
 925                                 success = minDiff;
 926                         }
 927
 928                         if (alignment != NULL) {  delete alignment;  }
 929
 930                 }
 931
 932
 933                 return success;
 934
 935         }
 936         catch(exception& e) {
 937                 m->errorOut(e, "TrimOligos", "stripLinker");
 938                 exit(1);
 939         }
 940 }
 941 //******************************************************************/
 942 bool TrimOligos::stripLinker(Sequence& seq){
 943         try {
 944
 945                 string rawSequence = seq.getUnaligned();
 946                 bool success = ldiffs +1;       //guilty until proven innocent
 947
 948                 for(int i=0;i<linker.size();i++){
 949                         string oligo = linker[i];
 950
 951                         if(rawSequence.length() < oligo.length()){
 952                                 success = ldiffs +10;
 953                                 break;
 954                         }
 955
 956                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 957                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 958                                 success = 0;
 959                                 break;
 960                         }
 961                 }
 962
 963         //if you found the linker or if you don't want to allow for diffs
 964                 if ((ldiffs == 0) || (success == 0)) { return success;  }
 965
 966                 else { //try aligning and see if you can find it
 967
 968                         int maxLength = 0;
 969
 970                         Alignment* alignment;
 971                         if (linker.size() > 0) {
 972                                 for(int i = 0; i < linker.size(); i++){
 973                                         if(linker[i].length() > maxLength){
 974                                                 maxLength = linker[i].length();
 975                                         }
 976                                 }
 977                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+ldiffs+1));
 978
 979                         }else{ alignment = NULL; }
 980
 981                         //can you find the barcode
 982                         int minDiff = 1e6;
 983                         int minCount = 1;
 984                         int minPos = 0;
 985
 986                         for(int i = 0; i < linker.size(); i++){
 987                                 string oligo = linker[i];
 988                                 //                              int length = oligo.length();
 989
 990                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 991                                         success = ldiffs + 10;
 992                                         break;
 993                                 }
 994
 995                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 996                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
 997                                 oligo = alignment->getSeqAAln();
 998                                 string temp = alignment->getSeqBAln();
 999
1000                                 int alnLength = oligo.length();
1001
1002                                 for(int i=oligo.length()-1;i>=0;i--){
1003                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
1004                                 }
1005                                 oligo = oligo.substr(0,alnLength);
1006                                 temp = temp.substr(0,alnLength);
1007
1008                                 int numDiff = countDiffs(oligo, temp);
1009
1010                                 if(numDiff < minDiff){
1011                                         minDiff = numDiff;
1012                                         minCount = 1;
1013                                         minPos = 0;
1014                                         for(int i=0;i<alnLength;i++){
1015                                                 if(temp[i] != '-'){
1016                                                         minPos++;
1017                                                 }
1018                                         }
1019                                 }
1020                                 else if(numDiff == minDiff){
1021                                         minCount++;
1022                                 }
1023
1024                         }
1025
1026                         if(minDiff > ldiffs)    {       success = minDiff;              }       //no good matches
1027                         else if(minCount > 1)   {       success = ldiffs + 100; }       //can't tell the difference between multiple barcodes
1028                         else{                                                                                                   //use the best match
1029                                 seq.setUnaligned(rawSequence.substr(minPos));
1030                                 success = minDiff;
1031                         }
1032
1033                         if (alignment != NULL) {  delete alignment;  }
1034
1035                 }
1036
1037                 return success;
1038
1039         }
1040         catch(exception& e) {
1041                 m->errorOut(e, "TrimOligos", "stripLinker");
1042                 exit(1);
1043         }
1044 }
1045
1046 //******************************************************************/
1047 bool TrimOligos::stripSpacer(Sequence& seq, QualityScores& qual){
1048         try {
1049                 string rawSequence = seq.getUnaligned();
1050                 bool success = sdiffs+1;        //guilty until proven innocent
1051
1052                 for(int i=0;i<spacer.size();i++){
1053                         string oligo = spacer[i];
1054
1055                         if(rawSequence.length() < oligo.length()){
1056                                 success = sdiffs+10;
1057                                 break;
1058                         }
1059
1060                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
1061                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
1062                                 if(qual.getName() != ""){
1063                                         qual.trimQScores(oligo.length(), -1);
1064                                 }
1065                                 success = 0;
1066                                 break;
1067                         }
1068                 }
1069
1070         //if you found the spacer or if you don't want to allow for diffs
1071                 if ((sdiffs == 0) || (success == 0)) { return success;  }
1072
1073                 else { //try aligning and see if you can find it
1074
1075                         int maxLength = 0;
1076
1077                         Alignment* alignment;
1078                         if (spacer.size() > 0) {
1079                                 for(int i = 0; i < spacer.size(); i++){
1080                                         if(spacer[i].length() > maxLength){
1081                                                 maxLength = spacer[i].length();
1082                                         }
1083                                 }
1084                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+sdiffs+1));
1085
1086                         }else{ alignment = NULL; }
1087
1088                         //can you find the barcode
1089                         int minDiff = 1e6;
1090                         int minCount = 1;
1091                         int minPos = 0;
1092
1093                         for(int i = 0; i < spacer.size(); i++){
1094                                 string oligo = spacer[i];
1095                                 //                              int length = oligo.length();
1096
1097                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
1098                                         success = sdiffs + 10;
1099                                         break;
1100                                 }
1101
1102                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1103                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
1104                                 oligo = alignment->getSeqAAln();
1105                                 string temp = alignment->getSeqBAln();
1106
1107                                 int alnLength = oligo.length();
1108
1109                                 for(int i=oligo.length()-1;i>=0;i--){
1110                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
1111                                 }
1112                                 oligo = oligo.substr(0,alnLength);
1113                                 temp = temp.substr(0,alnLength);
1114
1115                                 int numDiff = countDiffs(oligo, temp);
1116
1117                                 if(numDiff < minDiff){
1118                                         minDiff = numDiff;
1119                                         minCount = 1;
1120                                         minPos = 0;
1121                                         for(int i=0;i<alnLength;i++){
1122                                                 if(temp[i] != '-'){
1123                                                         minPos++;
1124                                                 }
1125                                         }
1126                                 }
1127                                 else if(numDiff == minDiff){
1128                                         minCount++;
1129                                 }
1130
1131                         }
1132
1133                         if(minDiff > sdiffs)    {       success = minDiff;              }       //no good matches
1134                         else if(minCount > 1)   {       success = sdiffs + 100; }       //can't tell the difference between multiple barcodes
1135                         else{                                                                                                   //use the best match
1136                                 seq.setUnaligned(rawSequence.substr(minPos));
1137
1138                                 if(qual.getName() != ""){
1139                                         qual.trimQScores(minPos, -1);
1140                                 }
1141                                 success = minDiff;
1142                         }
1143
1144                         if (alignment != NULL) {  delete alignment;  }
1145
1146                 }
1147
1148
1149                 return success;
1150
1151         }
1152         catch(exception& e) {
1153                 m->errorOut(e, "TrimOligos", "stripSpacer");
1154                 exit(1);
1155         }
1156 }
1157 //******************************************************************/
1158 bool TrimOligos::stripSpacer(Sequence& seq){
1159         try {
1160
1161                 string rawSequence = seq.getUnaligned();
1162                 bool success = sdiffs+1;        //guilty until proven innocent
1163
1164                 for(int i=0;i<spacer.size();i++){
1165                         string oligo = spacer[i];
1166
1167                         if(rawSequence.length() < oligo.length()){
1168                                 success = sdiffs+10;
1169                                 break;
1170                         }
1171
1172                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
1173                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
1174                                 success = 0;
1175                                 break;
1176                         }
1177                 }
1178
1179         //if you found the spacer or if you don't want to allow for diffs
1180                 if ((sdiffs == 0) || (success == 0)) { return success;  }
1181
1182                 else { //try aligning and see if you can find it
1183
1184                         int maxLength = 0;
1185
1186                         Alignment* alignment;
1187                         if (spacer.size() > 0) {
1188                                 for(int i = 0; i < spacer.size(); i++){
1189                                         if(spacer[i].length() > maxLength){
1190                                                 maxLength = spacer[i].length();
1191                                         }
1192                                 }
1193                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+sdiffs+1));
1194
1195                         }else{ alignment = NULL; }
1196
1197                         //can you find the barcode
1198                         int minDiff = 1e6;
1199                         int minCount = 1;
1200                         int minPos = 0;
1201
1202                         for(int i = 0; i < spacer.size(); i++){
1203                                 string oligo = spacer[i];
1204                                 //                              int length = oligo.length();
1205
1206                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
1207                                         success = sdiffs + 10;
1208                                         break;
1209                                 }
1210
1211                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1212                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
1213                                 oligo = alignment->getSeqAAln();
1214                                 string temp = alignment->getSeqBAln();
1215
1216                                 int alnLength = oligo.length();
1217
1218                                 for(int i=oligo.length()-1;i>=0;i--){
1219                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
1220                                 }
1221                                 oligo = oligo.substr(0,alnLength);
1222                                 temp = temp.substr(0,alnLength);
1223
1224                                 int numDiff = countDiffs(oligo, temp);
1225
1226                                 if(numDiff < minDiff){
1227                                         minDiff = numDiff;
1228                                         minCount = 1;
1229                                         minPos = 0;
1230                                         for(int i=0;i<alnLength;i++){
1231                                                 if(temp[i] != '-'){
1232                                                         minPos++;
1233                                                 }
1234                                         }
1235                                 }
1236                                 else if(numDiff == minDiff){
1237                                         minCount++;
1238                                 }
1239
1240                         }
1241
1242                         if(minDiff > sdiffs)    {       success = minDiff;              }       //no good matches
1243                         else if(minCount > 1)   {       success = sdiffs + 100; }       //can't tell the difference between multiple barcodes
1244                         else{                                                                                                   //use the best match
1245                                 seq.setUnaligned(rawSequence.substr(minPos));
1246                                 success = minDiff;
1247                         }
1248
1249                         if (alignment != NULL) {  delete alignment;  }
1250
1251                 }
1252
1253                 return success;
1254
1255         }
1256         catch(exception& e) {
1257                 m->errorOut(e, "TrimOligos", "stripSpacer");
1258                 exit(1);
1259         }
1260 }
1261
1262 //******************************************************************/
1263 bool TrimOligos::compareDNASeq(string oligo, string seq){
1264         try {
1265                 bool success = 1;
1266                 int length = oligo.length();
1267
1268                 for(int i=0;i<length;i++){
1269
1270                         if(oligo[i] != seq[i]){
1271                                 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C')    {       success = 0;    }
1272                                 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                                {       success = 0;    }
1273                                 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                                    {       success = 0;    }
1274                                 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                                    {       success = 0;    }
1275                                 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                                    {       success = 0;    }
1276                                 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                                    {       success = 0;    }
1277                                 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                                    {       success = 0;    }
1278                                 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                                    {       success = 0;    }
1279                                 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
1280                                 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
1281                                 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))   {       success = 0;    }
1282                                 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))   {       success = 0;    }
1283
1284                                 if(success == 0)        {       break;   }
1285                         }
1286                         else{
1287                                 success = 1;
1288                         }
1289                 }
1290
1291                 return success;
1292         }
1293         catch(exception& e) {
1294                 m->errorOut(e, "TrimOligos", "compareDNASeq");
1295                 exit(1);
1296         }
1297
1298 }
1299 //********************************************************************/
1300 int TrimOligos::countDiffs(string oligo, string seq){
1301         try {
1302
1303                 int length = oligo.length();
1304                 int countDiffs = 0;
1305
1306                 for(int i=0;i<length;i++){
1307
1308                         if(oligo[i] != seq[i]){
1309                                 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C' || oligo[i] == '-' || oligo[i] == '.')      {       countDiffs++;   }
1310                                 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                                {       countDiffs++;   }
1311                                 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                                    {       countDiffs++;   }
1312                                 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                                    {       countDiffs++;   }
1313                                 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                                    {       countDiffs++;   }
1314                                 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                                    {       countDiffs++;   }
1315                                 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                                    {       countDiffs++;   }
1316                                 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                                    {       countDiffs++;   }
1317                                 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))   {       countDiffs++;   }
1318                                 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))   {       countDiffs++;   }
1319                                 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))   {       countDiffs++;   }
1320                                 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))   {       countDiffs++;   }
1321                         }
1322
1323                 }
1324
1325                 return countDiffs;
1326         }
1327         catch(exception& e) {
1328                 m->errorOut(e, "TrimOligos", "countDiffs");
1329                 exit(1);
1330         }
1331 }
1332 /********************************************************************/
1333
1334
1335