ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/GAlnExtend.cpp
(Generate patch)
# Line 784 | Line 784
784    delete ed_script_fwd;
785    return alninfo;
786   }
787 +
788 + GXAlnInfo* match_RightEnd(const char* seqa, int seqa_len, const char* seqb, int seqb_len,
789 +            SGreedyAlignMem* gxmem, int match_reward, int mismatch_penalty, int Xdrop) {
790 +  bool editscript=false;
791 +  #ifdef GDEBUG
792 +   editscript=true;
793 +   GMessage("==========> matching Right (3') end ========\n");
794 +  #endif
795 +  GList<GXSeed> rseeds(true,true,false);
796 +  GXBandSet* alnbands=collectSeeds_R(rseeds, seqa, seqa_len, seqb, seqb_len);
797 +  GList<GXAlnInfo> galns(true,true,false);
798 +  int max_top_bands=5;
799 +  int top_band_count=0;
800 +  GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
801 +  for (int b=0;b<alnbands->Count();b++) {
802 +     if (alnbands->Get(b)->score<4) break;
803 +     top_band_count++;
804 +     GXBand& band=*(alnbands->Get(b));
805 +     band.seeds.setSorted(cmpSeedScore);
806 +     anchor_seeds.Add(band.seeds.First());
807 +     if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
808 +     }
809 +  #ifdef GDEBUG
810 +  GMessage("::: Collected %d anchor seeds.\n",anchor_seeds.Count());
811 +  #endif
812 +
813 +  for (int i=0;i<anchor_seeds.Count();i++) {
814 +    GXSeed& aseed=*anchor_seeds[i];
815 +    int a1=aseed.a_ofs+(aseed.len>>1)+1;
816 +    int a2=aseed.b_ofs+(aseed.len>>1)+1;
817 +    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
818 +                            seqb, a2, seqb_len, galn_TrimRight, editscript, gxmem,
819 +                            match_reward, mismatch_penalty, Xdrop);
820 +    if (alninfo && alninfo->pid>73 && alninfo->sr>seqb_len-4 && alninfo->ql<4)
821 +             galns.AddIfNew(alninfo, true);
822 +        else delete alninfo;
823 +    }
824 +
825 +  #ifdef GDEBUG
826 +  for (int i=0;i<galns.Count();i++) {
827 +    GXAlnInfo* alninfo=galns[i];
828 +    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
829 +                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
830 +    if (alninfo->gapinfo!=NULL) {
831 +      GMessage("Alignment:\n");
832 +      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
833 +      }
834 +    }
835 +  #endif
836 +  //---- done
837 +  delete alnbands;
838 +  if (galns.Count()) {
839 +    GXAlnInfo* bestaln=galns.Shift();
840 +    return bestaln;
841 +    }
842 +  else return NULL;
843 + }
844 +
845 + GXAlnInfo* match_LeftEnd(const char* seqa, int seqa_len, const char* seqb, int seqb_len,
846 +          SGreedyAlignMem* gxmem, int match_reward, int mismatch_penalty, int Xdrop) {
847 +  bool editscript=false;
848 +  #ifdef GDEBUG
849 +   editscript=true;
850 +   GMessage("==========> matching Left (5') end ========\n");
851 +  #endif
852 +  GList<GXSeed> rseeds(true,true,false);
853 +  GXBandSet* alnbands=collectSeeds_L(rseeds, seqa, seqa_len, seqb, seqb_len);
854 +  GList<GXAlnInfo> galns(true,true,false);
855 +  int max_top_bands=5;
856 +  int top_band_count=0;
857 +  GList<GXSeed> anchor_seeds(cmpSeedDiag, NULL, true); //stores unique seeds per diagonal
858 +  for (int b=0;b<alnbands->Count();b++) {
859 +     if (alnbands->Get(b)->score<4) break;
860 +     top_band_count++;
861 +     GXBand& band=*(alnbands->Get(b));
862 +     band.seeds.setSorted(cmpSeedScore);
863 +     anchor_seeds.Add(band.seeds.First());
864 +     if (anchor_seeds.Count()>2 || top_band_count>max_top_bands) break;
865 +     }
866 +  #ifdef GDEBUG
867 +  GMessage("::: Collected %d anchor seeds.\n",anchor_seeds.Count());
868 +  #endif
869 +
870 +  for (int i=0;i<anchor_seeds.Count();i++) {
871 +    GXSeed& aseed=*anchor_seeds[i];
872 +    int a1=aseed.a_ofs+(aseed.len>>1)+1;
873 +    int a2=aseed.b_ofs+(aseed.len>>1)+1;
874 +    GXAlnInfo* alninfo=GreedyAlignRegion(seqa, a1, seqa_len,
875 +                            seqb, a2, seqb_len, galn_TrimLeft, editscript, gxmem,
876 +                            match_reward, mismatch_penalty, Xdrop-2);
877 +    if (alninfo && alninfo->pid>82 && alninfo->sl<3 && alninfo->qr>seqa_len-4)
878 +         //TODO: at 5' end we should use a higher pid threshold
879 +         // we could also decrease the X-drop value!
880 +            galns.AddIfNew(alninfo, true);
881 +       else delete alninfo;
882 +    }
883 +
884 +  #ifdef GDEBUG
885 +  for (int i=0;i<galns.Count();i++) {
886 +    GXAlnInfo* alninfo=galns[i];
887 +    GMessage("a(%d..%d) align to b(%d..%d), score=%d, pid=%4.2f\n", alninfo->ql, alninfo->qr,
888 +                         alninfo->sl, alninfo->sr, alninfo->score, alninfo->pid);
889 +    if (alninfo->gapinfo!=NULL) {
890 +      GMessage("Alignment:\n");
891 +      alninfo->gapinfo->printAlignment(stderr, seqa, seqa_len, seqb, seqb_len);
892 +      }
893 +    }
894 +  #endif
895 +  //---- done
896 +  delete alnbands;
897 +  if (galns.Count()) {
898 +    GXAlnInfo* bestaln=galns.Shift();
899 +    return bestaln;
900 +    }
901 +  else return NULL;
902 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines