ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GAlnExtend.cpp
(Generate patch)
# Line 604 | Line 604
604   const int a_dropoff_score=7;
605   const int a_min_score=12; //at least 6 bases full match
606  
607 < // ------------------ adapter matching - simple k-mer seed & extend, no indels for now
607 > // ------------------ adaptor matching - simple k-mer seed & extend, no indels for now
608   //when a k-mer match is found, simply try to extend the alignment using a drop-off scheme
609   //check minimum score and
610 < //for 3' adapter trimming:
610 > //for 3' adaptor trimming:
611   //     require that the right end of the alignment for either the adaptor OR the read must be
612   //     < 3 distance from its right end
613 < // for 5' adapter trimming:
613 > // for 5' adaptor trimming:
614   //     require that the left end of the alignment for either the adaptor OR the read must
615   //     be at coordinate < 3 from start
616  
# Line 683 | Line 683
683    int mmovh5=(a_l<b_l)? a_l : b_l;
684    if (maxscore>=a_min_score && mmovh3<2 && mmovh5<2) {
685       if (a_l<a_ovh3) {
686 <        //adapter closer to the left end (typical for 5' adapter)
686 >        //adaptor closer to the left end (typical for 5' adaptor)
687          l5=a_r+1;
688          l3=alen-1;
689          }
690        else {
691 <        //adapter matching at the right end (typical for 3' adapter)
691 >        //adaptor matching at the right end (typical for 3' adaptor)
692          l5=0;
693          l3=a_l-1;
694          }
# Line 700 | Line 700
700    return false;
701   }
702  
703 <
704 < GXBandSet* collectSeeds_R(GList<GXSeed>& seeds, const char* seqa, int a_len,
705 <        GVec<uint16>* amers[], const char* seqb, int b_len) {
706 < int bimin=GMAX(0,(b_len-a_len-6));
707 < GXSeedTable gx(a_len, b_len);
708 < GXBandSet* diagstrips=new GXBandSet(a_len, b_len); //set of overlapping 3-diagonal strips
709 < for (int bi=0;bi<=b_len-6;bi++) {
703 > /*
704 > GXBandSet* collectSeeds_R(GList<GXSeed>& seeds, GXSeqData& sd) {
705 > int bimin=GMAX(0,(sd.blen-sd.alen-6));
706 > GXSeedTable gx(sd.alen, sd.blen);
707 > GXBandSet* diagstrips=new GXBandSet(sd.alen, sd.blen); //set of overlapping 3-diagonal strips
708 > //for (int bi=0;bi<=b_len-6;bi++) {
709 > for (int bi=sd.blen-6;bi>=0;bi--) {
710     //for each 6-mer of seqb
711 <   uint16 bv = get6mer((char*)&seqb[bi]);
712 <   GVec<uint16>* alocs = amers[bv];
711 >   uint16 bv = get6mer((char*)&(sd.bseq[bi]));
712 >   GVec<uint16>* alocs = sd.amers[bv];
713     if (alocs==NULL) continue;
714     //extend each hit
715     for (int h=0;h<alocs->Count();h++) {
# Line 720 | Line 720
720             continue;
721           for (int i=0;i<6;i++)
722              gx.x(ai+i,bi+i)=1;
723 <         //see if we can extend to the right
724 <         int aix=ai+6;
725 <         int bix=bi+6;
723 >         //see if we can extend to the left
724 >         int aix=ai-1;
725 >         int bix=bi-1;
726           int len=6;
727 <         while (bix<b_len && aix<a_len && seqa[aix]==seqb[bix]) {
727 >         while (bix>=0 && aix>=0 && sd.aseq[aix]==sd.bseq[bix]) {
728                   gx.x(aix,bix)=1;
729 <                 aix++;bix++;
729 >                 ai=aix;
730 >                 bi=bix;
731 >                 aix--;bix--;
732                   len++;
733                   }
734 <        /*if (len>22) {
735 <                //heuristics: very likely the best we can get
734 >        if (len>sd.amlen) {
735 >                //heuristics: likely the best we can get
736                  //quick match shortcut
737                  diagstrips->qmatch=new GXSeed(ai,bi,len);
738                  return diagstrips;
739 <                } */
739 >                }
740           if (bi<bimin && len<9) continue; //skip middle seeds that are not high scoring enough
741           GXSeed* newseed=new GXSeed(ai,bi,len);
742           seeds.Add(newseed);
743           diagstrips->addSeed(newseed);//add it to all 3 adjacent diagonals
744           //special last resort terminal match to be used if no better alignment is there
745 <         if (ai<2 && bi+len>b_len-2 &&
745 >         if (ai<2 && bi+len>sd.blen-2 &&
746                   (!diagstrips->tmatch || diagstrips->tmatch->len<len)) diagstrips->tmatch=newseed;
747           }
748    }
# Line 750 | Line 752
752   diagstrips->setSorted(true); //sort by score
753   return diagstrips;
754   }
755 <
756 < GXBandSet* collectSeeds_L(GList<GXSeed>& seeds, const char* seqa, int a_len,
757 <        GVec<uint16>* amers[], const char* seqb, int b_len) {
758 < //overlap of left (5') end of seqb
757 < //hash the last 12 bases of seqa:
758 < int bimax=GMIN((a_len+2), (b_len-6));
755 > */
756 > GXBandSet* collectSeeds(GList<GXSeed>& seeds, GXSeqData& sd) {
757 > int bimax=GMIN((sd.alen+2), (sd.blen-6));
758 > int bimin=GMAX(0,(sd.blen-sd.alen-6)); //from collectSeeds_R
759   //gx.init(a_maxlen, b_maxlen);
760 < GXSeedTable gx(a_len, b_len);
761 < GXBandSet* diagstrips=new GXBandSet(a_len, b_len); //set of overlapping 3-diagonal strips
762 < for (int bi=0;bi<=b_len-6;bi++) {
760 > GXSeedTable gx(sd.alen, sd.blen);
761 > GXBandSet* diagstrips=new GXBandSet(sd.alen, sd.blen); //set of overlapping 3-diagonal strips
762 > for (int bi=0;bi<=sd.blen-6;bi++) {
763     //for each 6-mer of seqb
764 <   uint16 bv = get6mer((char*)&seqb[bi]);
765 <   GVec<uint16>* alocs = amers[bv];
764 >   uint16 bv = get6mer((char*) & (sd.bseq[bi]));
765 >   GVec<uint16>* alocs = sd.amers[bv];
766     if (alocs==NULL) continue;
767     //extend each hit
768     for (int h=0;h<alocs->Count();h++) {
# Line 777 | Line 777
777           int aix=ai+6;
778           int bix=bi+6;
779           int len=6;
780 <         while (bix<b_len && aix<a_len && seqa[aix]==seqb[bix]) {
780 >         while (bix<sd.blen && aix<sd.alen && sd.aseq[aix]==sd.bseq[bix]) {
781                   gx.x(aix,bix)=1;
782                   aix++;bix++;
783                   len++;
784                   }
785 <         /* if (len>22) {
785 >         if (len>sd.amlen) {
786                  //heuristics: very likely the best we can get
787                  //quick match shortcut
788                  diagstrips->qmatch=new GXSeed(ai,bi,len);
789                  return diagstrips;
790 <                } */
791 <         if (bi>bimax && len<9) continue; //skip middle seeds that are not high scoring enough
790 >                }
791 >         if (bi>bimax && bi<bimin && len<9)
792 >                 //skip mid-sequence seeds that are not high scoring
793 >             continue;
794 >
795           GXSeed* newseed=new GXSeed(ai,bi,len);
796           seeds.Add(newseed);
797           diagstrips->addSeed(newseed);//add it to all 3 adjacent diagonals
798       //special last resort terminal match to be used if no better alignment is there
799 <     if (bi<2 && ai+len>=a_len-1 &&
800 <                 (!diagstrips->tmatch || diagstrips->tmatch->len<len))
801 <                   diagstrips->tmatch=newseed;
799 >     if (bi<2 && ai+len>=sd.alen-1 &&
800 >                 (!diagstrips->tmatch_l || diagstrips->tmatch_l->len<len))
801 >                      diagstrips->tmatch_l=newseed;
802 >     //collectSeeds_R:
803 >         if (ai<2 && bi+len>sd.blen-2 &&
804 >                 (!diagstrips->tmatch_r || diagstrips->tmatch_r->len<len))
805 >                  diagstrips->tmatch_r=newseed;
806       }
807   } //for each 6-mer of the read
808   for (int i=0;i<diagstrips->Count();i++) {
# Line 865 | Line 872
872    const char* s=s_seq+s_alnstart;
873    int s_avail=s_max-s_alnstart;
874    if (penalty<0) penalty=-penalty;
875 <  int MIN_GREEDY_SCORE=5*reward; //minimum score for an alignment to be reported for 0 diffs
875 >  int MIN_GREEDY_SCORE=6*reward; //minimum score for an alignment to be reported for 0 diffs
876    GXAlnInfo* alninfo=NULL;
877    bool freeAlnMem=(gxmem==NULL);
878    if (freeAlnMem) {
# Line 875 | Line 882
882    int retscore = 0;
883    int numdiffs = 0;
884    if (trim!=NULL && trim->type==galn_TrimLeft) {
885 +    //intent: trimming the left side
886      if (editscript)
887         ed_script_rev=new GXEditScript();
888  
889      int numdiffs_l = GXGreedyExtend(s_seq, s_alnstart, q_seq, q_alnstart, true, xdrop,
890                     reward, penalty, s_ext_l, q_ext_l, *gxmem, ed_script_rev);
891      //check this extension here and bail out if it's not a good extension
892 <    if (s_alnstart+1-s_ext_l>trim->boundary) {
892 >    if (s_ext_l+(trim->seedlen>>1) < trim->safelen &&
893 >        q_alnstart+1-q_ext_l>1 &&
894 >        s_alnstart+1-s_ext_l>trim->l_boundary) {
895        delete ed_script_rev;
896        if (freeAlnMem) delete gxmem;
897        return NULL;
# Line 906 | Line 916
916      //check extension here and bail out if not a good right extension
917      //assuming s_max is really at the right end of s_seq
918      if (trim!=NULL && trim->type==galn_TrimRight &&
919 <           s_alnstart+s_ext_r<trim->boundary) {
919 >        s_ext_r+(trim->seedlen>>1) < trim->safelen &&
920 >            q_alnstart+q_ext_r<q_max-2 &&
921 >            s_alnstart+s_ext_r<trim->r_boundary) {
922        delete ed_script_fwd;
923        if (freeAlnMem) delete gxmem;
924        return NULL;
# Line 961 | Line 973
973       reward, penalty, xdrop, gxmem, trim, editscript);
974   }
975  
976 < GXAlnInfo* match_RightEnd(const char* seqa, int seqa_len, GVec<uint16>* amers[],
977 <                 const char* seqb, int seqb_len, CGreedyAlignData* gxmem, int min_pid) {
976 > GXAlnInfo* match_adaptor(GXSeqData& sd, GAlnTrimType trim_type,
977 >                                     CGreedyAlignData* gxmem, int min_pid) {
978    bool editscript=false;
979    #ifdef GDEBUG
980     editscript=true;
981 <   GMessage("==========> matching Right (3') end : %s\n", seqa);
981 >   if (trim_type==galn_TrimLeft) {
982 >         GMessage("=======> searching left (5') end : %s\n", sd.aseq);
983 >     }
984 >   else if (trim_type==galn_TrimRight) {
985 >     GMessage("=======> searching right(3') end : %s\n", sd.aseq);
986 >     }
987 >   else if (trim_type==galn_TrimEither) {
988 >     GMessage("==========> searching  both ends : %s\n", sd.aseq);
989 >     }
990    #endif
991 <
972 <  CAlnTrim trimInfo(galn_TrimRight, seqb, seqb_len);
991 >  CAlnTrim trimInfo(trim_type, sd.bseq, sd.blen, sd.alen, sd.amlen);
992    GList<GXSeed> rseeds(true,true,false);
993 <  GXBandSet* alnbands=collectSeeds_R(rseeds, seqa, seqa_len, amers, seqb, seqb_len);
993 >  GXBandSet* alnbands=collectSeeds(rseeds, sd);
994    GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
995    //did we find a shortcut?
996    if (alnbands->qmatch) {
# Line 993 | Line 1012
1012         GXBand& band=*(alnbands->Get(b));
1013         band.seeds.setSorted(cmpSeedScore);
1014         anchor_seeds.Add(band.seeds.First());
1015 <       band.tested=true;
1015 >       //band.tested=true;
1016         if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
1017         }
1018      //#ifdef GDEBUG
# Line 1005 | Line 1024
1024      GXSeed& aseed=*anchor_seeds[i];
1025      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1026      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1027 <    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
1028 <                            seqb, a2, seqb_len, gxmem, &trimInfo, editscript);
1029 <    if (alninfo && alninfo->pid>=min_pid &&
1030 <        trimInfo.validate(seqa_len, alninfo->sl, alninfo->sr, alninfo->pid, alninfo->ql-1))
1027 >    trimInfo.seedlen=aseed.len;
1028 > #ifdef GDEBUG
1029 >    GMessage("\t::: align from seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1030 >           aseed.len);
1031 > #endif
1032 >    GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1033 >                            sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1034 >    if (alninfo && alninfo->pid>=min_pid && trimInfo.validate(alninfo))
1035               galns.AddIfNew(alninfo, true);
1036          else delete alninfo;
1037      }
1015  if (galns.Count()==0 && alnbands->tmatch) {
1016      //last resort seed
1017      GXSeed& aseed=*alnbands->tmatch;
1018      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1019      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1020      GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
1021                                seqb, a2, seqb_len, gxmem, &trimInfo, editscript);
1022      if (alninfo && alninfo->pid>=min_pid &&
1023           trimInfo.validate(seqa_len, alninfo->sl, alninfo->sr, alninfo->pid, alninfo->ql-1))
1024                 galns.AddIfNew(alninfo, true);
1025            else delete alninfo;
1026      }
1038  
1028  /*
1029  //special 3' end case: due to the seed scoring scheme being biased towards the 5' end of the read,
1030  //we should also try some seeds closer to the 3' end
1039    if (galns.Count()==0) {
1040 <    anchor_seeds.Clear();
1041 <    alnbands->setSorted(cmpDiagBands_R);
1042 <    int max_top_bands=4;
1043 <    int top_band_count=0;
1044 <    //#ifdef GDEBUG
1045 <    //GMessage(":::>> Retrying adjusting sort order.\n");
1046 <    //#endif
1047 <    if (alnbands->tmatch) {
1048 <      //anchor_seeds.setSorted(false);
1049 <      anchor_seeds.Add(alnbands->tmatch);
1050 <      }
1051 <    for (int b=0;b<alnbands->Count();b++) {
1052 <       if (alnbands->Get(b)->score<4) break;
1053 <       //#ifdef GDEBUG
1054 <       //GMessage("\tBand %d score: %d\n", b, alnbands->Get(b)->score);
1055 <       //#endif
1056 <       if (alnbands->Get(b)->tested) continue;
1057 <       top_band_count++;
1058 <       GXBand& band=*(alnbands->Get(b));
1059 <       band.seeds.setSorted(cmpSeedScore);
1060 <       anchor_seeds.Add(band.seeds.First());
1061 <       if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
1062 <       }
1063 <    //#ifdef GDEBUG
1064 <    //GMessage("::: Collected %d anchor seeds.\n",anchor_seeds.Count());
1065 <    //#endif
1058 <    for (int i=0;i<anchor_seeds.Count();i++) {
1059 <      GXSeed& aseed=*anchor_seeds[i];
1060 <      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1061 <      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1062 <      GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
1063 <                              seqb, a2, seqb_len, gxmem, &trimInfo, editscript);
1064 <      if (alninfo && alninfo->pid>=min_pid && trimInfo.validate(alninfo->sl, alninfo->sr, alninfo->pid, alninfo->ql-1))
1065 <               galns.AddIfNew(alninfo, true);
1066 <          else delete alninfo;
1040 >         //last resort: look for weaker terminal seeds
1041 >          GPVec<GXSeed> tmatches(2,false);
1042 >          if (trim_type!=galn_TrimRight) {
1043 >                 if (alnbands->tmatch_l)
1044 >                    tmatches.Add(alnbands->tmatch_l);
1045 >             }
1046 >          if (trim_type!=galn_TrimLeft) {
1047 >                 if (alnbands->tmatch_r)
1048 >                    tmatches.Add(alnbands->tmatch_r);
1049 >             }
1050 >          for (int i=0;i<tmatches.Count();i++) {
1051 >                GXSeed& aseed=*tmatches[i];
1052 >                int halfseed=aseed.len>>1;
1053 >                int a1=aseed.a_ofs+halfseed+1;
1054 >                int a2=aseed.b_ofs+halfseed+1;
1055 >                trimInfo.seedlen=aseed.len;
1056 > #ifdef GDEBUG
1057 >    GMessage("\t::: align from terminal seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1058 >           aseed.len);
1059 > #endif
1060 >        GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1061 >                                sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1062 >        if (alninfo && alninfo->pid>=min_pid && trimInfo.validate(alninfo))
1063 >                 galns.AddIfNew(alninfo, true);
1064 >             else delete alninfo;
1065 >        }//for each terminal seed
1066        }
1067 <    } */
1069 <
1070 <  //---- done
1067 >  //---- found all alignments
1068    delete alnbands;
1069 +  /*
1070 +  #ifdef GDEBUG
1071 +  //print all valid alignments found
1072 +  for (int i=0;i<galns.Count();i++) {
1073 +    GXAlnInfo* alninfo=galns[i];
1074 +    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
1075 +                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
1076 +    if (alninfo->gapinfo!=NULL) {
1077 +      GMessage("Alignment:\n");
1078 +      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
1079 +      }
1080 +    }
1081 +  #endif
1082 +  */
1083    if (galns.Count()) {
1084      GXAlnInfo* bestaln=galns.Shift();
1085      #ifdef GDEBUG
1086        GMessage("Best alignment: a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", bestaln->ql, bestaln->qr,
1087            bestaln->sl, bestaln->sr, bestaln->score, bestaln->pid);
1088        if (bestaln->gapinfo!=NULL) {
1089 <        bestaln->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
1089 >        bestaln->gapinfo->printAlignment(stderr, sd.aseq, sd.alen, sd.bseq, sd.blen);
1090          }
1091      #endif
1081
1092      return bestaln;
1093      }
1094    else return NULL;
1095   }
1096 <
1097 < GXAlnInfo* match_LeftEnd(const char* seqa, int seqa_len, GVec<uint16>* amers[],
1088 <                 const char* seqb, int seqb_len, CGreedyAlignData* gxmem, int min_pid) {
1096 > /*
1097 > GXAlnInfo* match_Left(GXSeqData& sd, CGreedyAlignData* gxmem, int min_pid) {
1098    bool editscript=false;
1099    #ifdef GDEBUG
1100     editscript=true;
1101 <   GMessage("==========> matching Left (5') end : %s\n", seqa);
1101 >   GMessage("==========> matching Left (5') end : %s\n", sd.aseq);
1102    #endif
1103 <  CAlnTrim trimInfo(galn_TrimLeft, seqb, seqb_len);
1103 >  CAlnTrim trimInfo(galn_TrimLeft, sd.bseq, sd.blen, sd.alen, sd.amlen);
1104    GList<GXSeed> rseeds(true,true,false);
1105 <  GXBandSet* alnbands=collectSeeds_L(rseeds, seqa, seqa_len, amers, seqb, seqb_len);
1105 >  GXBandSet* alnbands = collectSeeds(rseeds, sd);
1106    GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
1107    if (alnbands->qmatch) {
1108      #ifdef GDEBUG
# Line 1114 | Line 1123
1123         GXBand& band=*(alnbands->Get(b));
1124         band.seeds.setSorted(cmpSeedScore);
1125         anchor_seeds.Add(band.seeds.First());
1126 <       band.tested=true;
1126 >       //band.tested=true;
1127         if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
1128         }
1129      //#ifdef GDEBUG
# Line 1126 | Line 1135
1135      GXSeed& aseed=*anchor_seeds[i];
1136      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1137      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1138 <    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
1139 <                            seqb, a2, seqb_len, gxmem, &trimInfo, editscript);
1140 <    if (alninfo && alninfo->pid>=min_pid
1141 <           && trimInfo.validate(seqa_len, alninfo->sl, alninfo->sr, alninfo->pid, seqa_len-alninfo->qr))
1138 >    trimInfo.seedlen=aseed.len;
1139 > #ifdef GDEBUG
1140 >    GMessage("\t::: align from seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1141 >           aseed.len);
1142 > #endif
1143 >    GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1144 >                            sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1145 >    if (alninfo && alninfo->pid>=min_pid && trimInfo.validate(alninfo))
1146              galns.AddIfNew(alninfo, true);
1147         else delete alninfo;
1148      }
1149 <  if (galns.Count()==0 && alnbands->tmatch) {
1149 >  if (galns.Count()==0 && alnbands->tmatch_l) {
1150        //last resort seed
1151 <      GXSeed& aseed=*alnbands->tmatch;
1151 >      GXSeed& aseed=*alnbands->tmatch_l;
1152        int a1=aseed.a_ofs+(aseed.len>>1)+1;
1153        int a2=aseed.b_ofs+(aseed.len>>1)+1;
1154 <      GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
1155 <                              seqb, a2, seqb_len, gxmem, &trimInfo, editscript);
1156 <      if (alninfo && alninfo->pid>=min_pid &&
1157 <        trimInfo.validate(seqa_len, alninfo->sl, alninfo->sr, alninfo->pid, seqa_len-alninfo->qr))
1154 >      trimInfo.seedlen=aseed.len;
1155 >      GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1156 >                              sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1157 >      if (alninfo && alninfo->pid>=min_pid && trimInfo.validate(alninfo))
1158           galns.Add(alninfo);
1159        }
1147  /*
1148  #ifdef GDEBUG
1149  //print valid alignments found
1150  for (int i=0;i<galns.Count();i++) {
1151    GXAlnInfo* alninfo=galns[i];
1152    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
1153                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
1154    if (alninfo->gapinfo!=NULL) {
1155      GMessage("Alignment:\n");
1156      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
1157      }
1158    }
1159  #endif
1160  */
1160    //---- done
1161    delete alnbands;
1162    if (galns.Count()) {
# Line 1166 | Line 1165
1165        GMessage("Best alignment: a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", bestaln->ql, bestaln->qr,
1166            bestaln->sl, bestaln->sr, bestaln->score, bestaln->pid);
1167        if (bestaln->gapinfo!=NULL) {
1168 <        bestaln->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
1168 >        bestaln->gapinfo->printAlignment(stderr, sd.aseq, sd.alen,
1169 >                                          sd.bseq, sd.blen);
1170          }
1171      #endif
1172      return bestaln;
1173      }
1174    else return NULL;
1175   }
1176 + */

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines