ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/fqtrim/fqtrim.cpp
(Generate patch)
# Line 957 | Line 957
957   int wlen=rlen;
958   GXSeqData seqdata;
959   int numruns=revCompl ? 2 : 1;
960 + GList<GXAlnInfo> bestalns(true, true, false);
961   for (int ai=0;ai<adaptors3.Count();ai++) {
962     for (int r=0;r<numruns;r++) {
963       if (r) {
# Line 967 | Line 968
968              seqdata.update(adaptors3[ai].seq.chars(), adaptors3[ai].seq.length(),
969                   adaptors3[ai].pz, wseq.chars(), wlen, adaptors3[ai].amlen);
970          }
971 <
972 <  GXAlnInfo* bestaln=match_adaptor(seqdata, adaptors3[ai].trim_type, gxmem_r, 86);
973 <  if (bestaln) {
974 <     trimmed=true;
975 <     //keep unmatched region on the left OR right (the longer one)
976 <     if (bestaln->sl > wlen-bestaln->sr) {
977 <         //keep left side
978 <         l3-=(wlen-bestaln->sl+1);
979 <         if (l3<0) l3=0;
979 <         }
980 <     else { //keep right side
981 <         l5+=bestaln->sr;
982 <         if (l5>=rlen) l5=rlen-1;
983 <         }
984 <     delete bestaln;
985 <     if (l3-l5+1<min_read_len) return true;
986 <     wseq=seq.substr(l5,l3-l5+1);
987 <     wlen=wseq.length();
988 <     }
971 >     GXAlnInfo* aln=match_adaptor(seqdata, adaptors3[ai].trim_type, gxmem_r, 86);
972 >         if (aln) {
973 >           if (aln->strong) {
974 >                   trimmed=true;
975 >                   bestalns.Add(aln);
976 >                   break; //will check the rest next time
977 >                   }
978 >            else bestalns.Add(aln);
979 >           }
980     }//forward and reverse adaptors
981 +   if (trimmed) break; //will check the rest in the next cycle
982    }//for each 3' adaptor
983 <  return trimmed;
983 > if (bestalns.Count()>0) {
984 >           GXAlnInfo* aln=bestalns[0];
985 >           if (aln->sl-1 > wlen-aln->sr) {
986 >                   //keep left side
987 >                   l3-=(wlen-aln->sl+1);
988 >                   if (l3<0) l3=0;
989 >                   }
990 >           else { //keep right side
991 >                   l5+=aln->sr;
992 >                   if (l5>=rlen) l5=rlen-1;
993 >                   }
994 >           //delete aln;
995 >           //if (l3-l5+1<min_read_len) return true;
996 >           wseq=seq.substr(l5,l3-l5+1);
997 >           wlen=wseq.length();
998 >           return true; //break the loops here to report a good find
999 >     }
1000 >  return false;
1001   }
1002  
1003   bool trim_adaptor5(GStr& seq, int&l5, int &l3) {
# Line 1001 | Line 1010
1010   int wlen=rlen;
1011   GXSeqData seqdata;
1012   int numruns=revCompl ? 2 : 1;
1013 + GList<GXAlnInfo> bestalns(true, true, false);
1014   for (int ai=0;ai<adaptors5.Count();ai++) {
1015     for (int r=0;r<numruns;r++) {
1016       if (r) {
# Line 1011 | Line 1021
1021              seqdata.update(adaptors5[ai].seq.chars(), adaptors5[ai].seq.length(),
1022                   adaptors5[ai].pz, wseq.chars(), wlen, adaptors5[ai].amlen);
1023          }
1024 <         GXAlnInfo* bestaln=match_adaptor(seqdata, adaptors5[ai].trim_type, gxmem_l, 90);
1025 <         if (bestaln) {
1026 <           trimmed=true;
1027 <           if (bestaln->sl-1 > wlen-bestaln->sr) {
1024 >         GXAlnInfo* aln=match_adaptor(seqdata, adaptors5[ai].trim_type, gxmem_l, 90);
1025 >         if (aln) {
1026 >           if (aln->strong) {
1027 >                   trimmed=true;
1028 >                   bestalns.Add(aln);
1029 >                   break; //will check the rest next time
1030 >                   }
1031 >            else bestalns.Add(aln);
1032 >           }
1033 >         } //forward and reverse?
1034 >   if (trimmed) break; //will check the rest in the next cycle
1035 >  }//for each 5' adaptor
1036 >  if (bestalns.Count()>0) {
1037 >           GXAlnInfo* aln=bestalns[0];
1038 >           if (aln->sl-1 > wlen-aln->sr) {
1039                     //keep left side
1040 <                   l3-=(wlen-bestaln->sl+1);
1040 >                   l3-=(wlen-aln->sl+1);
1041                     if (l3<0) l3=0;
1042                     }
1043             else { //keep right side
1044 <                   l5+=bestaln->sr;
1044 >                   l5+=aln->sr;
1045                     if (l5>=rlen) l5=rlen-1;
1046                     }
1047 <           delete bestaln;
1048 <           if (l3-l5+1<min_read_len) return true;
1047 >           //delete aln;
1048 >           //if (l3-l5+1<min_read_len) return true;
1049             wseq=seq.substr(l5,l3-l5+1);
1050             wlen=wseq.length();
1051 <           }
1052 <         } //forward and reverse?
1053 <  }//for each 5' adaptor
1033 <  return trimmed;
1051 >           return true; //break the loops here to report a good find
1052 >     }
1053 >  return false;
1054   }
1055  
1056   //convert qvs to/from phred64 from/to phread33

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines