ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GAlnExtend.cpp
(Generate patch)
# Line 700 | Line 700
700    return false;
701   }
702  
703 <
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);
# Line 752 | Line 752
752   diagstrips->setSorted(true); //sort by score
753   return diagstrips;
754   }
755 <
756 < GXBandSet* collectSeeds_L(GList<GXSeed>& seeds, GXSeqData& sd) {
757 < //overlap of left (5') end of seqb
758 < //hash the last 12 bases of seqa:
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(sd.alen, sd.blen);
761   GXBandSet* diagstrips=new GXBandSet(sd.alen, sd.blen); //set of overlapping 3-diagonal strips
# Line 790 | Line 789
789                  return diagstrips;
790                  }
791           if (bi>bimax && len<9) continue; //skip middle seeds that are not high scoring enough
792 +         if (bi<bimin && len<9) continue; //collectSeeds_R
793 +
794           GXSeed* newseed=new GXSeed(ai,bi,len);
795           seeds.Add(newseed);
796           diagstrips->addSeed(newseed);//add it to all 3 adjacent diagonals
797       //special last resort terminal match to be used if no better alignment is there
798       if (bi<2 && ai+len>=sd.alen-1 &&
799 <                 (!diagstrips->tmatch || diagstrips->tmatch->len<len))
800 <                   diagstrips->tmatch=newseed;
799 >                 (!diagstrips->tmatch_l || diagstrips->tmatch_l->len<len))
800 >                      diagstrips->tmatch_l=newseed;
801 >     //collectSeeds_R:
802 >         if (ai<2 && bi+len>sd.blen-2 &&
803 >                 (!diagstrips->tmatch_r || diagstrips->tmatch_r->len<len))
804 >                  diagstrips->tmatch_r=newseed;
805       }
806   } //for each 6-mer of the read
807   for (int i=0;i<diagstrips->Count();i++) {
# Line 866 | Line 871
871    const char* s=s_seq+s_alnstart;
872    int s_avail=s_max-s_alnstart;
873    if (penalty<0) penalty=-penalty;
874 <  int MIN_GREEDY_SCORE=5*reward; //minimum score for an alignment to be reported for 0 diffs
874 >  int MIN_GREEDY_SCORE=6*reward; //minimum score for an alignment to be reported for 0 diffs
875    GXAlnInfo* alninfo=NULL;
876    bool freeAlnMem=(gxmem==NULL);
877    if (freeAlnMem) {
# Line 883 | Line 888
888      int numdiffs_l = GXGreedyExtend(s_seq, s_alnstart, q_seq, q_alnstart, true, xdrop,
889                     reward, penalty, s_ext_l, q_ext_l, *gxmem, ed_script_rev);
890      //check this extension here and bail out if it's not a good extension
891 <    if (s_ext_l+(trim->seedlen>>1) < trim->safelen && s_alnstart+1-s_ext_l>trim->boundary) {
891 >    if (s_ext_l+(trim->seedlen>>1) < trim->safelen &&
892 >        q_alnstart+1-q_ext_l>trim->boundary &&
893 >        s_alnstart+1-s_ext_l>trim->boundary) {
894        delete ed_script_rev;
895        if (freeAlnMem) delete gxmem;
896        return NULL;
# Line 909 | Line 916
916      //assuming s_max is really at the right end of s_seq
917      if (trim!=NULL && trim->type==galn_TrimRight &&
918          s_ext_r+(trim->seedlen>>1) < trim->safelen &&
919 +            q_alnstart+q_ext_r<q_max-3 &&
920              s_alnstart+s_ext_r<trim->boundary) {
921        delete ed_script_fwd;
922        if (freeAlnMem) delete gxmem;
# Line 964 | Line 972
972       reward, penalty, xdrop, gxmem, trim, editscript);
973   }
974  
975 < GXAlnInfo* match_RightEnd(GXSeqData& sd, CGreedyAlignData* gxmem,
968 <                 int min_pid) {
975 > GXAlnInfo* match_RightEnd(GXSeqData& sd, CGreedyAlignData* gxmem, int min_pid) {
976    bool editscript=false;
977    #ifdef GDEBUG
978     editscript=true;
# Line 974 | Line 981
981    CAlnTrim trimInfo(galn_TrimRight, sd.bseq, sd.blen, sd.amlen);
982    GList<GXSeed> rseeds(true,true,false);
983  
984 <  GXBandSet* alnbands=collectSeeds_R(rseeds, sd);
984 >  GXBandSet* alnbands=collectSeeds(rseeds, sd);
985    GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
986    //did we find a shortcut?
987    if (alnbands->qmatch) {
# Line 1009 | Line 1016
1016      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1017      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1018      trimInfo.seedlen=aseed.len;
1019 + #ifdef GDEBUG
1020 +    GMessage("\t::: align from seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1021 +           aseed.len);
1022 + #endif
1023      GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1024                              sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1025      if (alninfo && alninfo->pid>=min_pid &&
# Line 1016 | Line 1027
1027               galns.AddIfNew(alninfo, true);
1028          else delete alninfo;
1029      }
1030 <  if (galns.Count()==0 && alnbands->tmatch) {
1030 >  if (galns.Count()==0 && alnbands->tmatch_r) {
1031        //last resort seed
1032 <      GXSeed& aseed=*alnbands->tmatch;
1032 >      GXSeed& aseed=*alnbands->tmatch_r;
1033        int halfseed=aseed.len>>1;
1034        int a1=aseed.a_ofs+halfseed+1;
1035        int a2=aseed.b_ofs+halfseed+1;
1036        trimInfo.seedlen=aseed.len;
1037 + #ifdef GDEBUG
1038 +    GMessage("\t::: align from seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1039 +           aseed.len);
1040 + #endif
1041        GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1042                                  sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1043        if (alninfo && alninfo->pid>=min_pid &&
# Line 1098 | Line 1113
1113    #endif
1114    CAlnTrim trimInfo(galn_TrimLeft, sd.bseq, sd.blen, sd.amlen);
1115    GList<GXSeed> rseeds(true,true,false);
1116 <  GXBandSet* alnbands=collectSeeds_L(rseeds, sd);
1116 >  GXBandSet* alnbands = collectSeeds(rseeds, sd);
1117    GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
1118    if (alnbands->qmatch) {
1119      #ifdef GDEBUG
# Line 1132 | Line 1147
1147      int a1=aseed.a_ofs+(aseed.len>>1)+1;
1148      int a2=aseed.b_ofs+(aseed.len>>1)+1;
1149      trimInfo.seedlen=aseed.len;
1150 + #ifdef GDEBUG
1151 +    GMessage("\t::: align from seed (%d, %d)of len %d.\n",aseed.a_ofs, aseed.b_ofs,
1152 +           aseed.len);
1153 + #endif
1154      GXAlnInfo* alninfo=GreedyAlignRegion(sd.aseq, a1, sd.alen,
1155                              sd.bseq, a2, sd.blen, gxmem, &trimInfo, editscript);
1156      if (alninfo && alninfo->pid>=min_pid
# Line 1140 | Line 1159
1159              galns.AddIfNew(alninfo, true);
1160         else delete alninfo;
1161      }
1162 <  if (galns.Count()==0 && alnbands->tmatch) {
1162 >  if (galns.Count()==0 && alnbands->tmatch_l) {
1163        //last resort seed
1164 <      GXSeed& aseed=*alnbands->tmatch;
1164 >      GXSeed& aseed=*alnbands->tmatch_l;
1165        int a1=aseed.a_ofs+(aseed.len>>1)+1;
1166        int a2=aseed.b_ofs+(aseed.len>>1)+1;
1167        trimInfo.seedlen=aseed.len;

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines