开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
sequence.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：137k
源码类别：

生物技术

开发平台：

C/C++

sequence.cpp：源码内容
  1.     case eOverlap_CheckIntervals:
  2.     case eOverlap_Interval:
  3.         revert_locations = true;
  4.         // there's no break here - proceed to "default"
  5.     default:
  6.         // Require intervals overlap
  7.         annot_overlap_type = SAnnotSelector::eOverlap_Intervals;
  8.         break;
  9.     }
  11.     CConstRef<CSeq_feat> feat_ref;
  12.     int diff = -1;
  14.     CFeat_CI feat_it(scope, loc, SAnnotSelector()
  15.         .SetFeatType(feat_type)
  16.         .SetFeatSubtype(feat_subtype)
  17.         .SetOverlapType(annot_overlap_type)
  18.         .SetResolveTSE());
  19.     for ( ; feat_it; ++feat_it) {
  20.         // treat subset as a special case
  21.         int cur_diff = ( !revert_locations ) ?
  22.             TestForOverlap(loc, feat_it->GetLocation(), overlap_type) :
  23.             TestForOverlap(feat_it->GetLocation(), loc, overlap_type);
  24.         if (cur_diff < 0)
  25.             continue;
  26.         if ( cur_diff < diff  ||  diff < 0 ) {
  27.             diff = cur_diff;
  28.             feat_ref = &feat_it->GetMappedFeature();
  29.         }
  30.     }
  31.     return feat_ref;
  32. }
  35. CConstRef<CSeq_feat> GetBestOverlappingFeat(const CSeq_loc& loc,
  36.                                             CSeqFeatData::E_Choice feat_type,
  37.                                             EOverlapType overlap_type,
  38.                                             CScope& scope)
  39. {
  40.     return x_GetBestOverlappingFeat(loc,
  41.         feat_type, CSeqFeatData::eSubtype_any,
  42.         overlap_type, scope);
  43. }
  46. CConstRef<CSeq_feat> GetBestOverlappingFeat(const CSeq_loc& loc,
  47.                                             CSeqFeatData::ESubtype feat_type,
  48.                                             EOverlapType overlap_type,
  49.                                             CScope& scope)
  50. {
  51.     return x_GetBestOverlappingFeat(loc,
  52.         CSeqFeatData::GetTypeFromSubtype(feat_type), feat_type,
  53.         overlap_type, scope);
  54. }
  57. CConstRef<CSeq_feat> GetOverlappingGene(const CSeq_loc& loc, CScope& scope)
  58. {
  59.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_gene,
  60.         eOverlap_Contains, scope);
  61. }
  64. CConstRef<CSeq_feat> GetOverlappingmRNA(const CSeq_loc& loc, CScope& scope)
  65. {
  66.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_mRNA,
  67.         eOverlap_Contains, scope);
  68. }
  71. CConstRef<CSeq_feat> GetOverlappingCDS(const CSeq_loc& loc, CScope& scope)
  72. {
  73.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_cdregion,
  74.         eOverlap_Contains, scope);
  75. }
  78. CConstRef<CSeq_feat> GetOverlappingPub(const CSeq_loc& loc, CScope& scope)
  79. {
  80.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_pub,
  81.         eOverlap_Contains, scope);
  82. }
  85. CConstRef<CSeq_feat> GetOverlappingSource(const CSeq_loc& loc, CScope& scope)
  86. {
  87.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_biosrc,
  88.         eOverlap_Contains, scope);
  89. }
  92. CConstRef<CSeq_feat> GetOverlappingOperon(const CSeq_loc& loc, CScope& scope)
  93. {
  94.     return GetBestOverlappingFeat(loc, CSeqFeatData::eSubtype_operon,
  95.         eOverlap_Contains, scope);
  96. }
  99. int SeqLocPartialCheck(const CSeq_loc& loc, CScope* scope)
  100. {
  101.     unsigned int retval = 0;
  102.     if (!scope) {
  103.         return retval;
  104.     }
  106.     // Find first and last Seq-loc
  107.     const CSeq_loc *first = 0, *last = 0;
  108.     for ( CSeq_loc_CI loc_iter(loc); loc_iter; ++loc_iter ) {
  109.         if ( first == 0 ) {
  110.             first = &(loc_iter.GetSeq_loc());
  111.         }
  112.         last = &(loc_iter.GetSeq_loc());
  113.     }
  114.     if (!first) {
  115.         return retval;
  116.     }
  118.     CSeq_loc_CI i2(loc, CSeq_loc_CI::eEmpty_Allow);
  119.     for ( ; i2; ++i2 ) {
  120.         const CSeq_loc* slp = &(i2.GetSeq_loc());
  121.         switch (slp->Which()) {
  122.         case CSeq_loc::e_Null:
  123.             if (slp == first) {
  124.                 retval |= eSeqlocPartial_Start;
  125.             } else if (slp == last) {
  126.                 retval |= eSeqlocPartial_Stop;
  127.             } else {
  128.                 retval |= eSeqlocPartial_Internal;
  129.             }
  130.             break;
  131.         case CSeq_loc::e_Int:
  132.             {
  133.                 const CSeq_interval& itv = slp->GetInt();
  134.                 if (itv.IsSetFuzz_from()) {
  135.                     const CInt_fuzz& fuzz = itv.GetFuzz_from();
  136.                     if (fuzz.Which() == CInt_fuzz::e_Lim) {
  137.                         CInt_fuzz::ELim lim = fuzz.GetLim();
  138.                         if (lim == CInt_fuzz::eLim_gt) {
  139.                             retval |= eSeqlocPartial_Limwrong;
  140.                         } else if (lim == CInt_fuzz::eLim_lt  ||
  141.                             lim == CInt_fuzz::eLim_unk) {
  142.                             if (itv.IsSetStrand()  &&
  143.                                 itv.GetStrand() == eNa_strand_minus) {
  144.                                 if (slp == last) {
  145.                                     retval |= eSeqlocPartial_Stop;
  146.                                 } else {
  147.                                     retval |= eSeqlocPartial_Internal;
  148.                                 }
  149.                                 if (itv.GetFrom() != 0) {
  150.                                     if (slp == last) {
  151.                                         retval |= eSeqlocPartial_Nostop;
  152.                                     } else {
  153.                                         retval |= eSeqlocPartial_Nointernal;
  154.                                     }
  155.                                 }
  156.                             } else {
  157.                                 if (slp == first) {
  158.                                     retval |= eSeqlocPartial_Start;
  159.                                 } else {
  160.                                     retval |= eSeqlocPartial_Internal;
  161.                                 }
  162.                                 if (itv.GetFrom() != 0) {
  163.                                     if (slp == first) {
  164.                                         retval |= eSeqlocPartial_Nostart;
  165.                                     } else {
  166.                                         retval |= eSeqlocPartial_Nointernal;
  167.                                     }
  168.                                 }
  169.                             }
  170.                         }
  171.                     }
  172.                 }
  173.                 
  174.                 if (itv.IsSetFuzz_to()) {
  175.                     const CInt_fuzz& fuzz = itv.GetFuzz_to();
  176.                     CInt_fuzz::ELim lim = fuzz.IsLim() ? 
  177.                         fuzz.GetLim() : CInt_fuzz::eLim_unk;
  178.                     if (lim == CInt_fuzz::eLim_lt) {
  179.                         retval |= eSeqlocPartial_Limwrong;
  180.                     } else if (lim == CInt_fuzz::eLim_gt  ||
  181.                         lim == CInt_fuzz::eLim_unk) {
  182.                         CBioseq_Handle hnd =
  183.                             scope->GetBioseqHandle(itv.GetId());
  184.                         bool miss_end = false;
  185.                         if ( hnd ) {
  186.                             CBioseq_Handle::TBioseqCore bc = hnd.GetBioseqCore();
  187.                             
  188.                             if (itv.GetTo() != bc->GetInst().GetLength() - 1) {
  189.                                 miss_end = true;
  190.                             }
  191.                         }
  192.                         if (itv.IsSetStrand()  &&
  193.                             itv.GetStrand() == eNa_strand_minus) {
  194.                             if (slp == first) {
  195.                                 retval |= eSeqlocPartial_Start;
  196.                             } else {
  197.                                 retval |= eSeqlocPartial_Internal;
  198.                             }
  199.                             if (miss_end) {
  200.                                 if (slp == first /* was last */) {
  201.                                     retval |= eSeqlocPartial_Nostart;
  202.                                 } else {
  203.                                     retval |= eSeqlocPartial_Nointernal;
  204.                                 }
  205.                             }
  206.                         } else {
  207.                             if (slp == last) {
  208.                                 retval |= eSeqlocPartial_Stop;
  209.                             } else {
  210.                                 retval |= eSeqlocPartial_Internal;
  211.                             }
  212.                             if (miss_end) {
  213.                                 if (slp == last) {
  214.                                     retval |= eSeqlocPartial_Nostop;
  215.                                 } else {
  216.                                     retval |= eSeqlocPartial_Nointernal;
  217.                                 }
  218.                             }
  219.                         }
  220.                     }
  221.                 }
  222.             }
  223.             break;
  224.         case CSeq_loc::e_Pnt:
  225.             if (slp->GetPnt().IsSetFuzz()) {
  226.                 const CInt_fuzz& fuzz = slp->GetPnt().GetFuzz();
  227.                 if (fuzz.Which() == CInt_fuzz::e_Lim) {
  228.                     CInt_fuzz::ELim lim = fuzz.GetLim();
  229.                     if (lim == CInt_fuzz::eLim_gt  ||
  230.                         lim == CInt_fuzz::eLim_lt  ||
  231.                         lim == CInt_fuzz::eLim_unk) {
  232.                         if (slp == first) {
  233.                             retval |= eSeqlocPartial_Start;
  234.                         } else if (slp == last) {
  235.                             retval |= eSeqlocPartial_Stop;
  236.                         } else {
  237.                             retval |= eSeqlocPartial_Internal;
  238.                         }
  239.                     }
  240.                 }
  241.             }
  242.             break;
  243.         case CSeq_loc::e_Packed_pnt:
  244.             if (slp->GetPacked_pnt().IsSetFuzz()) {
  245.                 const CInt_fuzz& fuzz = slp->GetPacked_pnt().GetFuzz();
  246.                 if (fuzz.Which() == CInt_fuzz::e_Lim) {
  247.                     CInt_fuzz::ELim lim = fuzz.GetLim();
  248.                     if (lim == CInt_fuzz::eLim_gt  ||
  249.                         lim == CInt_fuzz::eLim_lt  ||
  250.                         lim == CInt_fuzz::eLim_unk) {
  251.                         if (slp == first) {
  252.                             retval |= eSeqlocPartial_Start;
  253.                         } else if (slp == last) {
  254.                             retval |= eSeqlocPartial_Stop;
  255.                         } else {
  256.                             retval |= eSeqlocPartial_Internal;
  257.                         }
  258.                     }
  259.                 }
  260.             }
  261.             break;
  262.         case CSeq_loc::e_Whole:
  263.         {
  264.             // Get the Bioseq referred to by Whole
  265.             CBioseq_Handle bsh = scope->GetBioseqHandle(slp->GetWhole());
  266.             if ( !bsh ) {
  267.                 break;
  268.             }
  269.             // Check for CMolInfo on the biodseq
  270.             CSeqdesc_CI di( bsh, CSeqdesc::e_Molinfo );
  271.             if ( !di ) {
  272.                 // If no CSeq_descr, nothing can be done
  273.                 break;
  274.             }
  275.             // First try to loop through CMolInfo
  276.             const CMolInfo& mi = di->GetMolinfo();
  277.             if (!mi.IsSetCompleteness()) {
  278.                 continue;
  279.             }
  280.             switch (mi.GetCompleteness()) {
  281.             case CMolInfo::eCompleteness_no_left:
  282.                 if (slp == first) {
  283.                     retval |= eSeqlocPartial_Start;
  284.                 } else {
  285.                     retval |= eSeqlocPartial_Internal;
  286.                 }
  287.                 break;
  288.             case CMolInfo::eCompleteness_no_right:
  289.                 if (slp == last) {
  290.                     retval |= eSeqlocPartial_Stop;
  291.                 } else {
  292.                     retval |= eSeqlocPartial_Internal;
  293.                 }
  294.                 break;
  295.             case CMolInfo::eCompleteness_partial:
  296.                 retval |= eSeqlocPartial_Other;
  297.                 break;
  298.             case CMolInfo::eCompleteness_no_ends:
  299.                 retval |= eSeqlocPartial_Start;
  300.                 retval |= eSeqlocPartial_Stop;
  301.                 break;
  302.             default:
  303.                 break;
  304.             }
  305.             break;
  306.         }
  307.         default:
  308.             break;
  309.         }
  310.     }
  311.     return retval;
  312. }
  315. typedef CSeq_loc::TRange    TRange;
  316. typedef vector<TRange>      TRangeVec;
  319. CSeq_loc* SeqLocMerge
  320. (const CBioseq_Handle& target,
  321.  const CSeq_loc& loc1,
  322.  const CSeq_loc& loc2,
  323.  TSeqLocFlags flags)
  324. {
  325.     _ASSERT(target);
  327.     vector< CConstRef<CSeq_loc> > locs;
  328.     locs.push_back(CConstRef<CSeq_loc>(&loc1));
  329.     locs.push_back(CConstRef<CSeq_loc>(&loc2));
  330.     return SeqLocMerge(target, locs, flags);
  331. }
  334. static void s_CollectRanges(const CSeq_loc& loc, TRangeVec& ranges)
  335. {
  336.     // skipping empty / nulls
  337.     for ( CSeq_loc_CI li(loc); li; ++li ) {
  338.         ranges.push_back(li.GetRange());
  339.     }
  340. }
  343. static void s_MergeRanges(TRangeVec& ranges, TSeqLocFlags flags)
  344. {
  345.     bool merge = (flags & fMergeIntervals) != 0;
  346.     bool fuse  = (flags & fFuseAbutting) != 0;
  348.     if ( !merge  &&  !fuse ) {
  349.         return;
  350.     }
  352.     TRangeVec::iterator next = ranges.begin();
  353.     TRangeVec::iterator curr = next++;
  354.     bool advance = true;
  355.     while ( next != ranges.end() ) {
  356.         advance = true;
  357.         if ( (merge  &&  curr->GetTo() >= next->GetFrom())  ||
  358.              (fuse   &&  curr->GetTo() + 1 == next->GetFrom()) ) {
  359.             *curr += *next;
  360.             next = ranges.erase(next);
  361.             advance = false;
  362.         }
  363.         if ( advance ) {
  364.             curr = next++;
  365.         }
  366.     }
  367. }
  370. static CSeq_loc* s_CreateSeqLoc
  371. (CSeq_id& id,
  372.  const TRangeVec& ranges,
  373.  ENa_strand strand,
  374.  bool rearranged,
  375.  bool add_null)
  376. {
  377.     static const  CSeq_loc null_loc(CSeq_loc::e_Null);
  379.     bool has_pnt = false;
  380.     bool has_int = false;
  382.     if ( ranges.size() < 2 ) {
  383.         add_null = false;
  384.     }
  386.     if ( !add_null ) {
  387.         ITERATE(TRangeVec, it, ranges) {
  388.             if ( it->GetLength() == 1 ) {
  389.                 has_pnt = true;
  390.             } else {
  391.                 has_int = true;
  392.             }
  393.         }
  394.     }
  395.     
  396.     CSeq_loc::E_Choice loc_type = CSeq_loc::e_not_set;
  397.     if ( add_null ) {
  398.         loc_type = CSeq_loc::e_Mix;
  399.     } else if ( has_pnt  &&  has_int ) {  // points and intervals  => mix location
  400.         loc_type = CSeq_loc::e_Mix;
  401.     } else if ( has_pnt ) {  // only points
  402.         loc_type = ranges.size() == 1 ? CSeq_loc::e_Pnt : CSeq_loc::e_Packed_pnt;
  403.     } else if ( has_int ) {  // only intervals
  404.         loc_type = ranges.size() == 1 ? CSeq_loc::e_Int : CSeq_loc::e_Packed_int;
  405.     } else {
  406.         loc_type = CSeq_loc::e_Null;
  407.     }
  409.     CRef<CSeq_loc> result;
  410.     switch ( loc_type ) {
  411.     case CSeq_loc::e_Pnt:
  412.         {
  413.             result.Reset(new CSeq_loc(id, ranges.front().GetFrom(), strand));
  414.             break;
  415.         }
  416.     case CSeq_loc::e_Packed_pnt:
  417.         {
  418.             CPacked_seqpnt::TPoints points(ranges.size());;
  419.             ITERATE(TRangeVec, it, ranges) {
  420.                 points.push_back(it->GetFrom());
  421.             }
  422.             result.Reset(new CSeq_loc(id, points, strand));
  423.             
  424.             break;
  425.         }
  426.     case CSeq_loc::e_Int:
  427.         {
  428.             const TRange& r = ranges.front();
  429.             result.Reset(new CSeq_loc(id, r.GetFrom(), r.GetTo(), strand));
  430.             break;
  431.         }
  432.     case CSeq_loc::e_Packed_int:
  433.         {
  434.             result.Reset(new CSeq_loc(id, ranges, strand));
  435.             if ( rearranged ) {
  436.                 // the first 2 intervals span the origin
  437.                 CPacked_seqint::Tdata ivals = result->SetPacked_int().Set();
  438.                 CPacked_seqint::Tdata::iterator it = ivals.begin();
  439.                 if ( strand == eNa_strand_minus ) {
  440.                     (*it)->SetFuzz_from().SetLim(CInt_fuzz::eLim_circle);
  441.                     ++it;
  442.                     (*it)->SetFuzz_to().SetLim(CInt_fuzz::eLim_circle);
  443.                 } else {
  444.                     (*it)->SetFuzz_to().SetLim(CInt_fuzz::eLim_circle);
  445.                     ++it;
  446.                     (*it)->SetFuzz_from().SetLim(CInt_fuzz::eLim_circle);
  447.                 }
  448.             }
  449.             break;
  450.         }
  451.     case CSeq_loc::e_Mix:
  452.         {
  453.             result.Reset(new CSeq_loc);
  454.             CSeq_loc_mix& mix = result->SetMix();
  455.             TRangeVec::const_iterator last = ranges.end(); --last;
  456.             ITERATE(TRangeVec, it, ranges) {
  457.                 if ( it->GetLength() == 1 ) {
  458.                     mix.AddSeqLoc(*new CSeq_loc(id, it->GetFrom(), strand));
  459.                 } else {
  460.                     mix.AddSeqLoc(*new CSeq_loc(id, it->GetFrom(), it->GetTo(), strand));
  461.                 }
  462.                 if ( add_null  &&  it != last ) {
  463.                     mix.AddSeqLoc(null_loc);
  464.                 }
  465.             }
  466.             break;
  467.         }
  468.     case CSeq_loc::e_Null:
  469.         {
  470.             result.Reset(new CSeq_loc);
  471.             result->SetNull();
  472.             break;
  473.         }
  474.     default:
  475.         {
  476.             result.Reset();
  477.         }
  478.     };
  480.     return result.Release();
  481. }
  484. static void s_RearrangeRanges(TRangeVec& ranges)
  485. {
  486.     deque<TRange> temp(ranges.size());
  487.     copy(ranges.begin(), ranges.end(), temp.begin());
  488.     ranges.clear();
  489.     temp.push_front(temp.back());  
  490.     temp.pop_back();
  491.     copy(temp.begin(), temp.end(), back_inserter(ranges));
  492. }
  495. CSeq_loc* SeqLocMergeOne
  496. (const CBioseq_Handle& target,
  497.  const CSeq_loc& loc,
  498.  TSeqLocFlags flags)
  499. {
  500.     _ASSERT(target);
  502.     TRangeVec ranges;
  504.     CRef<CSeq_id> id(new CSeq_id);
  505.     id->Assign(GetId(target, eGetId_Best));
  507.     const CSeq_inst& inst = target.GetInst();
  508.     CSeq_inst::TTopology topology = inst.CanGetTopology() ?
  509.         inst.GetTopology() : CSeq_inst::eTopology_not_set;
  510.     bool circular = (topology == CSeq_inst::eTopology_circular);
  511.     if ( circular ) {  // circular topology overrides fSingleInterval flag
  512.         flags &= ~fSingleInterval;
  513.     }
  514.     TSeqPos seq_len = inst.IsSetLength() ? inst.GetLength() : 0;
  516.     // map the location to the target bioseq
  517.     CRef<CSeq_loc> mapped_loc(target.MapLocation(loc));
  518.     _ASSERT(IsOneBioseq(*mapped_loc));  // doesn't have multiple bioseqs
  520.     ENa_strand strand = GetStrand(*mapped_loc);
  521.     bool rearranged = false;
  523.     if ( (flags & fSingleInterval) != 0 ) {
  524.             ranges.push_back(mapped_loc->GetTotalRange());
  525.     } else {
  526.         if ( strand == eNa_strand_other ) {  // multiple strands in location
  527.             return mapped_loc.Release();
  528.         }
  529.     
  530.         s_CollectRanges(*mapped_loc, ranges);
  531.         sort(ranges.begin(), ranges.end());  // on the plus strand
  532.         s_MergeRanges(ranges, flags);
  533.         if ( strand == eNa_strand_minus ) {
  534.             reverse(ranges.begin(), ranges.end());
  535.         }
  536.         // if circular bioseq and ranges span the origin, rearrange
  537.         if ( circular ) {
  538.             if ( ((strand == eNa_strand_minus)             &&
  539.                   (ranges.front().GetTo() == seq_len - 1)  &&
  540.                   (ranges.back().GetFrom() == 0))           ||
  541.                  ((strand != eNa_strand_minus)             &&
  542.                   (ranges.front().GetFrom() == 0)          &&
  543.                   (ranges.back().GetTo() == seq_len - 1)) ) {
  544.                 rearranged = true;
  545.                 s_RearrangeRanges(ranges);
  546.             }
  547.         }
  548.     }   
  550.     return s_CreateSeqLoc(*id, ranges, strand, rearranged, (flags & fAddNulls) != 0);
  551. }
  554. static CSeq_interval* s_SeqIntRevCmp(const CSeq_interval& i, CScope* scope)
  555. {
  556.     CRef<CSeq_interval> rev_int(new CSeq_interval);
  557.     rev_int->Assign(i);
  558.     
  559.     ENa_strand s = i.CanGetStrand() ? i.GetStrand() : eNa_strand_unknown;
  560.     rev_int->SetStrand(Reverse(s));
  562.     return rev_int.Release();
  563. }
  566. static CSeq_point* s_SeqPntRevCmp(const CSeq_point& pnt, CScope* scope)
  567. {
  568.     CRef<CSeq_point> rev_pnt(new CSeq_point);
  569.     rev_pnt->Assign(pnt);
  570.     
  571.     ENa_strand s = pnt.CanGetStrand() ? pnt.GetStrand() : eNa_strand_unknown;
  572.     rev_pnt->SetStrand(Reverse(s));
  574.     return rev_pnt.Release();
  575. }
  578. CSeq_loc* SeqLocRevCmp(const CSeq_loc& loc, CScope* scope)
  579. {
  580.     CRef<CSeq_loc> rev_loc(new CSeq_loc);
  582.     switch ( loc.Which() ) {
  584.     // -- reverse is the same.
  585.     case CSeq_loc::e_Null:
  586.     case CSeq_loc::e_Empty:
  587.     case CSeq_loc::e_Whole:
  588.         rev_loc->Assign(loc);
  589.         break;
  591.     // -- just reverse the strand
  592.     case CSeq_loc::e_Int:
  593.         rev_loc->SetInt(*s_SeqIntRevCmp(loc.GetInt(), scope));
  594.         break;
  595.     case CSeq_loc::e_Pnt:
  596.         rev_loc->SetPnt(*s_SeqPntRevCmp(loc.GetPnt(), scope));
  597.         break;
  598.     case CSeq_loc::e_Packed_pnt:
  599.         rev_loc->SetPacked_pnt().Assign(loc.GetPacked_pnt());
  600.         rev_loc->SetPacked_pnt().SetStrand(Reverse(GetStrand(loc, scope)));
  601.         break;
  603.     // -- possibly more than one sequence
  604.     case CSeq_loc::e_Packed_int:
  605.     {
  606.         // reverse each interval and store them in reverse order
  607.         typedef CRef< CSeq_interval > TInt;
  608.         CPacked_seqint& pint = rev_loc->SetPacked_int();
  609.         ITERATE (CPacked_seqint::Tdata, it, loc.GetPacked_int().Get()) {
  610.             pint.Set().push_front(TInt(s_SeqIntRevCmp(**it, scope)));
  611.         }
  612.         break;
  613.     }
  614.     case CSeq_loc::e_Mix:
  615.     {
  616.         // reverse each location and store them in reverse order
  617.         typedef CRef< CSeq_loc > TLoc;
  618.         CSeq_loc_mix& mix = rev_loc->SetMix();
  619.         ITERATE (CSeq_loc_mix::Tdata, it, loc.GetMix().Get()) {
  620.             mix.Set().push_front(TLoc(SeqLocRevCmp(**it, scope)));
  621.         }
  622.         break;
  623.     }
  624.     case CSeq_loc::e_Equiv:
  625.     {
  626.         // reverse each location (no need to reverse order)
  627.         typedef CRef< CSeq_loc > TLoc;
  628.         CSeq_loc_equiv& e = rev_loc->SetEquiv();
  629.         ITERATE (CSeq_loc_equiv::Tdata, it, loc.GetEquiv().Get()) {
  630.             e.Set().push_back(TLoc(SeqLocRevCmp(**it, scope)));
  631.         }
  632.         break;
  633.     }
  635.     case CSeq_loc::e_Bond:
  636.     {
  637.         CSeq_bond& bond = rev_loc->SetBond();
  638.         bond.SetA(*s_SeqPntRevCmp(loc.GetBond().GetA(), scope));
  639.         if ( loc.GetBond().CanGetB() ) {
  640.             bond.SetA(*s_SeqPntRevCmp(loc.GetBond().GetB(), scope));
  641.         }
  642.     }
  643.         
  644.     // -- not supported
  645.     case CSeq_loc::e_Feat:
  646.     default:
  647.         NCBI_THROW(CException, eUnknown,
  648.             "CSeq_loc::SeqLocRevCmp -- unsupported location type");
  649.     }
  651.     return rev_loc.Release();
  652. }
  655. // Get the encoding CDS feature of a given protein sequence.
  656. const CSeq_feat* GetCDSForProduct(const CBioseq& product, CScope* scope)
  657. {
  658.     if ( scope == 0 ) {
  659.         return 0;
  660.     }
  662.     return GetCDSForProduct(scope->GetBioseqHandle(product));
  663. }
  665. const CSeq_feat* GetCDSForProduct(const CBioseq_Handle& bsh)
  666. {
  667.     if ( bsh ) {
  668.         CFeat_CI fi(bsh, 
  669.                     0, 0,
  670.                     CSeqFeatData::e_Cdregion,
  671.                     SAnnotSelector::eOverlap_Intervals,
  672.                     SAnnotSelector::eResolve_TSE,
  673.                     CFeat_CI::e_Product,
  674.                     0);
  675.         if ( fi ) {
  676.             // return the first one (should be the one packaged on the
  677.             // nuc-prot set).
  678.             return &(fi->GetOriginalFeature());
  679.         }
  680.     }
  682.     return 0;
  683. }
  686. // Get the mature peptide feature of a protein
  687. const CSeq_feat* GetPROTForProduct(const CBioseq& product, CScope* scope)
  688. {
  689.     if ( scope == 0 ) {
  690.         return 0;
  691.     }
  693.     return GetPROTForProduct(scope->GetBioseqHandle(product));
  694. }
  696. const CSeq_feat* GetPROTForProduct(const CBioseq_Handle& bsh)
  697. {
  698.     if ( bsh ) {
  699.         CFeat_CI fi(bsh, 
  700.                     0, 0,
  701.                     CSeqFeatData::e_Prot,
  702.                     SAnnotSelector::eOverlap_Intervals,
  703.                     SAnnotSelector::eResolve_TSE,
  704.                     CFeat_CI::e_Product,
  705.                     0);
  706.         if ( fi ) {
  707.             return &(fi->GetOriginalFeature());
  708.         }
  709.     }
  711.     return 0;
  712. }
  716. // Get the encoding mRNA feature of a given mRNA (cDNA) bioseq.
  717. const CSeq_feat* GetmRNAForProduct(const CBioseq& product, CScope* scope)
  718. {
  719.     if ( scope == 0 ) {
  720.         return 0;
  721.     }
  723.     return GetmRNAForProduct(scope->GetBioseqHandle(product));
  724. }
  726. const CSeq_feat* GetmRNAForProduct(const CBioseq_Handle& bsh)
  727. {
  728.     if ( bsh ) {
  729.         SAnnotSelector as;
  730.         as.SetFeatSubtype(CSeqFeatData::eSubtype_mRNA);
  731.         as.SetByProduct();
  733.         CFeat_CI fi(bsh, 0, 0, as);
  734.         if ( fi ) {
  735.             return &(fi->GetOriginalFeature());
  736.         }
  737.     }
  739.     return 0;
  740. }
  743. // Get the encoding sequence of a protein
  744. const CBioseq* GetNucleotideParent(const CBioseq& product, CScope* scope)
  745. {
  746.     if ( scope == 0 ) {
  747.         return 0;
  748.     }
  749.     CBioseq_Handle bsh = GetNucleotideParent(scope->GetBioseqHandle(product));
  750.     return bsh ? &(bsh.GetBioseq()) : 0;
  751. }
  753. CBioseq_Handle GetNucleotideParent(const CBioseq_Handle& bsh)
  754. {
  755.     // If protein use CDS to get to the encoding Nucleotide.
  756.     // if nucleotide (cDNA) use mRNA feature.
  757.     const CSeq_feat* sfp = bsh.GetBioseqMolType() == CSeq_inst::eMol_aa ?
  758.         GetCDSForProduct(bsh) : GetmRNAForProduct(bsh);
  760.     CBioseq_Handle ret;
  761.     if ( sfp ) {
  762.         ret = bsh.GetScope().GetBioseqHandle(sfp->GetLocation());
  763.     }
  764.     return ret;
  765. }
  768. END_SCOPE(sequence)
  771. void CFastaOstream::Write(const CBioseq_Handle& handle,
  772.                           const CSeq_loc* location)
  773. {
  774.     WriteTitle(handle);
  775.     WriteSequence(handle, location);
  776. }
  779. void CFastaOstream::WriteTitle(const CBioseq_Handle& handle)
  780. {
  781.     m_Out << '>' << CSeq_id::GetStringDescr(*handle.GetBioseqCore(),
  782.                                             CSeq_id::eFormat_FastA)
  783.           << ' ' << sequence::GetTitle(handle) << NcbiEndl;
  784. }
  787. // XXX - replace with SRelLoc?
  788. struct SGap {
  789.     SGap(TSeqPos start, TSeqPos length) : m_Start(start), m_Length(length) { }
  790.     TSeqPos GetEnd(void) const { return m_Start + m_Length - 1; }
  792.     TSeqPos m_Start, m_Length;
  793. };
  794. typedef list<SGap> TGaps;
  797. static bool s_IsGap(const CSeq_loc& loc, CScope& scope)
  798. {
  799.     if (loc.IsNull()) {
  800.         return true;
  801.     }
  803.     CTypeConstIterator<CSeq_id> id(loc);
  804.     CBioseq_Handle handle = scope.GetBioseqHandle(*id);
  805.     if (handle  &&  handle.GetBioseqCore()->GetInst().GetRepr()
  806.         == CSeq_inst::eRepr_virtual) {
  807.         return true;
  808.     }
  810.     return false; // default
  811. }
  814. static TGaps s_FindGaps(const CSeq_ext& ext, CScope& scope)
  815. {
  816.     TSeqPos pos = 0;
  817.     TGaps   gaps;
  819.     switch (ext.Which()) {
  820.     case CSeq_ext::e_Seg:
  821.         ITERATE (CSeg_ext::Tdata, it, ext.GetSeg().Get()) {
  822.             TSeqPos length = sequence::GetLength(**it, &scope);
  823.             if (s_IsGap(**it, scope)) {
  824.                 gaps.push_back(SGap(pos, length));
  825.             }
  826.             pos += length;
  827.         }
  828.         break;
  830.     case CSeq_ext::e_Delta:
  831.         ITERATE (CDelta_ext::Tdata, it, ext.GetDelta().Get()) {
  832.             switch ((*it)->Which()) {
  833.             case CDelta_seq::e_Loc:
  834.             {
  835.                 const CSeq_loc& loc = (*it)->GetLoc();
  836.                 TSeqPos length = sequence::GetLength(loc, &scope);
  837.                 if (s_IsGap(loc, scope)) {
  838.                     gaps.push_back(SGap(pos, length));
  839.                 }
  840.                 pos += length;
  841.                 break;
  842.             }
  844.             case CDelta_seq::e_Literal:
  845.             {
  846.                 const CSeq_literal& lit    = (*it)->GetLiteral();
  847.                 TSeqPos             length = lit.GetLength();
  848.                 if ( !lit.IsSetSeq_data() ) {
  849.                     gaps.push_back(SGap(pos, length));
  850.                 }
  851.                 pos += length;
  852.                 break;
  853.             }
  855.             default:
  856.                 ERR_POST(Warning << "CFastaOstream::WriteSequence: "
  857.                          "unsupported Delta-seq selection "
  858.                          << CDelta_seq::SelectionName((*it)->Which()));
  859.                 break;
  860.             }
  861.         }
  863.     default:
  864.         break;
  865.     }
  867.     return gaps;
  868. }
  871. static TGaps s_AdjustGaps(const TGaps& gaps, const CSeq_loc& location)
  872. {
  873.     // assume location matches handle
  874.     const TSeqPos         kMaxPos = numeric_limits<TSeqPos>::max();
  875.     TSeqPos               pos     = 0;
  876.     TGaps::const_iterator gap_it  = gaps.begin();
  877.     TGaps                 adjusted_gaps;
  878.     SGap                  new_gap(kMaxPos, 0);
  880.     for (CSeq_loc_CI loc_it(location);  loc_it  &&  gap_it != gaps.end();
  881.          pos += loc_it.GetRange().GetLength(), ++loc_it) {
  882.         CSeq_loc_CI::TRange range = loc_it.GetRange();
  884.         if (new_gap.m_Start != kMaxPos) {
  885.             // in progress
  886.             if (gap_it->GetEnd() < range.GetFrom()) {
  887.                 adjusted_gaps.push_back(new_gap);
  888.                 new_gap.m_Start = kMaxPos;
  889.                 ++gap_it;
  890.             } else if (gap_it->GetEnd() <= range.GetTo()) {
  891.                 new_gap.m_Length += gap_it->GetEnd() - range.GetFrom() + 1;
  892.                 adjusted_gaps.push_back(new_gap);
  893.                 new_gap.m_Start = kMaxPos;
  894.                 ++gap_it;
  895.             } else {
  896.                 new_gap.m_Length += range.GetLength();
  897.                 continue;
  898.             }
  899.         }
  901.         while (gap_it->GetEnd() < range.GetFrom()) {
  902.             ++gap_it; // skip
  903.         }
  905.         if (gap_it->m_Start <= range.GetFrom()) {
  906.             if (gap_it->GetEnd() <= range.GetTo()) {
  907.                 adjusted_gaps.push_back
  908.                     (SGap(pos, gap_it->GetEnd() - range.GetFrom() + 1));
  909.                 ++gap_it;
  910.             } else {
  911.                 new_gap.m_Start  = pos;
  912.                 new_gap.m_Length = range.GetLength();
  913.                 continue;
  914.             }
  915.         }
  917.         while (gap_it->m_Start <= range.GetTo()) {
  918.             TSeqPos pos2 = pos + gap_it->m_Start - range.GetFrom();
  919.             if (gap_it->GetEnd() <= range.GetTo()) {
  920.                 adjusted_gaps.push_back(SGap(pos2, gap_it->m_Length));
  921.                 ++gap_it;
  922.             } else {
  923.                 new_gap.m_Start  = pos2;
  924.                 new_gap.m_Length = range.GetTo() - gap_it->m_Start + 1;
  925.             }
  926.         }
  927.     }
  929.     if (new_gap.m_Start != kMaxPos) {
  930.         adjusted_gaps.push_back(new_gap);
  931.     }
  933.     return adjusted_gaps;
  934. }
  937. void CFastaOstream::WriteSequence(const CBioseq_Handle& handle,
  938.                                   const CSeq_loc* location)
  939. {
  940.     CConstRef<CBioseq> seq  = handle.GetCompleteBioseq();
  941.     const CSeq_inst& inst = seq->GetInst();
  942.     if ( !(m_Flags & eAssembleParts)  &&  !inst.IsSetSeq_data() ) {
  943.         return;
  944.     }
  946.     CSeqVector v = (location
  947.                     ? handle.GetSequenceView(*location,
  948.                                              CBioseq_Handle::eViewMerged,
  949.                                              CBioseq_Handle::eCoding_Iupac)
  950.                     : handle.GetSeqVector(CBioseq_Handle::eCoding_Iupac));
  951.     bool is_na = inst.GetMol() != CSeq_inst::eMol_aa;
  952.     // autodetection is sometimes broken (!)
  953.     v.SetCoding(is_na ? CSeq_data::e_Iupacna : CSeq_data::e_Iupacaa);
  955.     TSeqPos              size = v.size();
  956.     TSeqPos              pos  = 0;
  957.     CSeqVector::TResidue gap  = v.GetGapChar();
  958.     string               buffer;
  959.     TGaps                gaps;
  960.     CScope&              scope   = handle.GetScope();
  961.     const TSeqPos        kMaxPos = numeric_limits<TSeqPos>::max();
  963.     if ( !inst.IsSetSeq_data()  &&  inst.IsSetExt() ) {
  964.         gaps = s_FindGaps(inst.GetExt(), scope);
  965.         if (location) {
  966.             gaps = s_AdjustGaps(gaps, *location);
  967.         }
  968.     }
  969.     gaps.push_back(SGap(kMaxPos, 0));
  971.     while (pos < size) {
  972.         unsigned int limit = min(m_Width,
  973.                                  min(size, gaps.front().m_Start) - pos);
  974.         v.GetSeqData(pos, pos + limit, buffer);
  975.         pos += limit;
  976.         if (limit > 0) {
  977.             m_Out << buffer << 'n';
  978.         }
  979.         if (pos == gaps.front().m_Start) {
  980.             if (m_Flags & eInstantiateGaps) {
  981.                 for (TSeqPos l = gaps.front().m_Length; l > 0; l -= m_Width) {
  982.                     m_Out << string(min(l, m_Width), gap) << 'n';
  983.                 }
  984.             } else {
  985.                 m_Out << "-n";
  986.             }
  987.             pos += gaps.front().m_Length;
  988.             gaps.pop_front();
  989.         }
  990.     }
  991.     m_Out << NcbiFlush;
  992. }
  995. void CFastaOstream::Write(CSeq_entry& entry, const CSeq_loc* location)
  996. {
  997.     CObjectManager om;
  998.     CScope         scope(om);
  1000.     scope.AddTopLevelSeqEntry(entry);
  1001.     for (CTypeConstIterator<CBioseq> it(entry);  it;  ++it) {
  1002.         if ( !SkipBioseq(*it) ) {
  1003.             CBioseq_Handle h = scope.GetBioseqHandle(*it->GetId().front());
  1004.             Write(h, location);
  1005.         }
  1006.     }
  1007. }
  1010. void CFastaOstream::Write(CBioseq& seq, const CSeq_loc* location)
  1011. {
  1012.     CSeq_entry entry;
  1013.     entry.SetSeq(seq);
  1014.     Write(entry, location);
  1015. }
  1019. /////////////////////////////////////////////////////////////////////////////
  1020. //
  1021. // sequence translation
  1022. //
  1025. template <class Container>
  1026. void x_Translate(const Container& seq,
  1027.                  string& prot,
  1028.                  const CGenetic_code* code)
  1029. {
  1030.     // reserve our space
  1031.     const size_t mod = seq.size() % 3;
  1032.     prot.erase();
  1033.     prot.reserve(seq.size() / 3 + (mod ? 1 : 0));
  1035.     // get appropriate translation table
  1036.     const CTrans_table & tbl =
  1037.         (code ? CGen_code_table::GetTransTable(*code) :
  1038.                 CGen_code_table::GetTransTable(1));
  1040.     // main loop through bases
  1041.     typename Container::const_iterator start = seq.begin();
  1043.     size_t i;
  1044.     size_t k;
  1045.     size_t state = 0;
  1046.     size_t length = seq.size() / 3;
  1047.     for (i = 0;  i < length;  ++i) {
  1049.         // loop through one codon at a time
  1050.         for (k = 0;  k < 3;  ++k, ++start) {
  1051.             state = tbl.NextCodonState(state, *start);
  1052.         }
  1054.         // save translated amino acid
  1055.         prot.append(1, tbl.GetCodonResidue(state));
  1056.     }
  1058.     if (mod) {
  1059.         LOG_POST(Warning <<
  1060.                  "translation of sequence whose length "
  1061.                  "is not an even number of codons");
  1062.         for (k = 0;  k < mod;  ++k, ++start) {
  1063.             state = tbl.NextCodonState(state, *start);
  1064.         }
  1066.         for ( ;  k < 3;  ++k) {
  1067.             state = tbl.NextCodonState(state, 'N');
  1068.         }
  1070.         // save translated amino acid
  1071.         prot.append(1, tbl.GetCodonResidue(state));
  1072.     }
  1073. }
  1076. void CSeqTranslator::Translate(const string& seq, string& prot,
  1077.                                const CGenetic_code* code,
  1078.                                bool include_stop,
  1079.                                bool remove_trailing_X)
  1080. {
  1081.     x_Translate(seq, prot, code);
  1082. }
  1085. void CSeqTranslator::Translate(const CSeqVector& seq, string& prot,
  1086.                                const CGenetic_code* code,
  1087.                                bool include_stop,
  1088.                                bool remove_trailing_X)
  1089. {
  1090.     x_Translate(seq, prot, code);
  1091. }
  1094. void CSeqTranslator::Translate(const CSeq_loc& loc,
  1095.                                const CBioseq_Handle& handle,
  1096.                                string& prot,
  1097.                                const CGenetic_code* code,
  1098.                                bool include_stop,
  1099.                                bool remove_trailing_X)
  1100. {
  1101.     CSeqVector seq =
  1102.         handle.GetSequenceView(loc, CBioseq_Handle::eViewConstructed,
  1103.                                CBioseq_Handle::eCoding_Iupac);
  1104.     x_Translate(seq, prot, code);
  1105. }
  1110. void CCdregion_translate::ReadSequenceByLocation (string& seq,
  1111.                                                   const CBioseq_Handle& bsh,
  1112.                                                   const CSeq_loc& loc)
  1114. {
  1115.     // get vector of sequence under location
  1116.     CSeqVector seqv = bsh.GetSequenceView (loc,
  1117.                                            CBioseq_Handle::eViewConstructed,
  1118.                                            CBioseq_Handle::eCoding_Iupac);
  1119.     seqv.GetSeqData(0, seqv.size(), seq);
  1120. }
  1122. void CCdregion_translate::TranslateCdregion (string& prot,
  1123.                                              const CBioseq_Handle& bsh,
  1124.                                              const CSeq_loc& loc,
  1125.                                              const CCdregion& cdr,
  1126.                                              bool include_stop,
  1127.                                              bool remove_trailing_X,
  1128.                                              bool* alt_start)
  1129. {
  1130.     // clear contents of result string
  1131.     prot.erase();
  1132.     if ( alt_start != 0 ) {
  1133.         *alt_start = false;
  1134.     }
  1136.     // copy bases from coding region location
  1137.     string bases = "";
  1138.     ReadSequenceByLocation (bases, bsh, loc);
  1140.     // calculate offset from frame parameter
  1141.     int offset = 0;
  1142.     if (cdr.IsSetFrame ()) {
  1143.         switch (cdr.GetFrame ()) {
  1144.             case CCdregion::eFrame_two :
  1145.                 offset = 1;
  1146.                 break;
  1147.             case CCdregion::eFrame_three :
  1148.                 offset = 2;
  1149.                 break;
  1150.             default :
  1151.                 break;
  1152.         }
  1153.     }
  1155.     int dnalen = bases.size () - offset;
  1156.     if (dnalen < 1) return;
  1158.     // pad bases string if last codon is incomplete
  1159.     bool incomplete_last_codon = false;
  1160.     int mod = dnalen % 3;
  1161.     if ( mod != 0 ) {
  1162.         incomplete_last_codon = true;
  1163.         bases += (mod == 1) ? "NN" : "N";
  1164.         dnalen += 3 - mod;
  1165.     }
  1166.     _ASSERT((dnalen >= 3)  &&  ((dnalen % 3) == 0));
  1168.     // resize output protein translation string
  1169.     prot.resize(dnalen / 3);
  1171.     // get appropriate translation table
  1172.     const CTrans_table& tbl = cdr.IsSetCode() ?
  1173.         CGen_code_table::GetTransTable(cdr.GetCode()) :
  1174.         CGen_code_table::GetTransTable(1);
  1176.     // main loop through bases
  1177.     string::const_iterator it = bases.begin();
  1178.     string::const_iterator end = bases.end();
  1179.     for ( int i = 0; i < offset; ++i ) {
  1180.         ++it;
  1181.     }
  1182.     unsigned int state = 0, j = 0;
  1183.     while ( it != end ) {
  1184.         // get one codon at a time
  1185.         state = tbl.NextCodonState(state, *it++);
  1186.         state = tbl.NextCodonState(state, *it++);
  1187.         state = tbl.NextCodonState(state, *it++);
  1189.         // save translated amino acid
  1190.         prot[j++] = tbl.GetCodonResidue(state);
  1191.     }
  1193.     // check for alternative start codon
  1194.     if ( alt_start != 0 ) {
  1195.         state = 0;
  1196.         state = tbl.NextCodonState (state, bases[0]);
  1197.         state = tbl.NextCodonState (state, bases[1]);
  1198.         state = tbl.NextCodonState (state, bases[2]);
  1199.         if ( tbl.IsAltStart(state) ) {
  1200.             *alt_start = true;
  1201.         }
  1202.     }
  1204.     // if complete at 5' end, require valid start codon
  1205.     if (offset == 0  &&  (! loc.IsPartialLeft ())) {
  1206.         state = tbl.SetCodonState (bases [offset], bases [offset + 1], bases [offset + 2]);
  1207.         prot [0] = tbl.GetStartResidue (state);
  1208.     }
  1210.     // code break substitution
  1211.     if (cdr.IsSetCode_break ()) {
  1212.         SIZE_TYPE protlen = prot.size ();
  1213.         ITERATE (CCdregion::TCode_break, code_break, cdr.GetCode_break ()) {
  1214.             const CRef <CCode_break> brk = *code_break;
  1215.             const CSeq_loc& cbk_loc = brk->GetLoc ();
  1216.             TSeqPos seq_pos = sequence::LocationOffset (loc, cbk_loc, sequence::eOffset_FromStart, &bsh.GetScope ());
  1217.             seq_pos -= offset;
  1218.             SIZE_TYPE i = seq_pos / 3;
  1219.             if (i >= 0 && i < protlen) {
  1220.               const CCode_break::C_Aa& c_aa = brk->GetAa ();
  1221.               if (c_aa.IsNcbieaa ()) {
  1222.                 prot [i] = c_aa.GetNcbieaa ();
  1223.               }
  1224.             }
  1225.         }
  1226.     }
  1228.     // optionally truncate at first terminator
  1229.     if (! include_stop) {
  1230.         SIZE_TYPE protlen = prot.size ();
  1231.         for (SIZE_TYPE i = 0; i < protlen; i++) {
  1232.             if (prot [i] == '*') {
  1233.                 prot.resize (i);
  1234.                 return;
  1235.             }
  1236.         }
  1237.     }
  1239.     // if padding was needed, trim ambiguous last residue
  1240.     if (incomplete_last_codon) {
  1241.         int protlen = prot.size ();
  1242.         if (protlen > 0 && prot [protlen - 1] == 'X') {
  1243.             protlen--;
  1244.             prot.resize (protlen);
  1245.         }
  1246.     }
  1248.     // optionally remove trailing X on 3' partial coding region
  1249.     if (remove_trailing_X) {
  1250.         int protlen = prot.size ();
  1251.         while (protlen > 0 && prot [protlen - 1] == 'X') {
  1252.             protlen--;
  1253.         }
  1254.         prot.resize (protlen);
  1255.     }
  1256. }
  1259. void CCdregion_translate::TranslateCdregion
  1260. (string& prot,
  1261.  const CSeq_feat& cds,
  1262.  CScope& scope,
  1263.  bool include_stop,
  1264.  bool remove_trailing_X,
  1265.  bool* alt_start)
  1266. {
  1267.     _ASSERT(cds.GetData().IsCdregion());
  1269.     prot.erase();
  1271.     CBioseq_Handle bsh = scope.GetBioseqHandle(cds.GetLocation());
  1272.     if ( !bsh ) {
  1273.         return;
  1274.     }
  1276.     CCdregion_translate::TranslateCdregion(
  1277.         prot,
  1278.         bsh,
  1279.         cds.GetLocation(),
  1280.         cds.GetData().GetCdregion(),
  1281.         include_stop,
  1282.         remove_trailing_X,
  1283.         alt_start);
  1284. }
  1287. SRelLoc::SRelLoc(const CSeq_loc& parent, const CSeq_loc& child, CScope* scope,
  1288.                  SRelLoc::TFlags flags)
  1289.     : m_ParentLoc(&parent)
  1290. {
  1291.     typedef CSeq_loc::TRange TRange0;
  1292.     for (CSeq_loc_CI cit(child);  cit;  ++cit) {
  1293.         const CSeq_id& cseqid  = cit.GetSeq_id();
  1294.         TRange0        crange  = cit.GetRange();
  1295.         if (crange.IsWholeTo()  &&  scope) {
  1296.             // determine actual end
  1297.             crange.SetToOpen(sequence::GetLength(cit.GetSeq_id(), scope));
  1298.         }
  1299.         ENa_strand     cstrand = cit.GetStrand();
  1300.         TSeqPos        pos     = 0;
  1301.         for (CSeq_loc_CI pit(parent);  pit;  ++pit) {
  1302.             ENa_strand pstrand = pit.GetStrand();
  1303.             TRange0    prange  = pit.GetRange();
  1304.             if (prange.IsWholeTo()  &&  scope) {
  1305.                 // determine actual end
  1306.                 prange.SetToOpen(sequence::GetLength(pit.GetSeq_id(), scope));
  1307.             }
  1308.             if ( !sequence::IsSameBioseq(cseqid, pit.GetSeq_id(), scope) ) {
  1309.                 pos += prange.GetLength();
  1310.                 continue;
  1311.             }
  1312.             CRef<TRange>         intersection(new TRange);
  1313.             TSeqPos              abs_from, abs_to;
  1314.             CConstRef<CInt_fuzz> fuzz_from, fuzz_to;
  1315.             if (crange.GetFrom() >= prange.GetFrom()) {
  1316.                 abs_from  = crange.GetFrom();
  1317.                 fuzz_from = cit.GetFuzzFrom();
  1318.                 if (abs_from == prange.GetFrom()) {
  1319.                     // subtract out parent fuzz, if any
  1320.                     const CInt_fuzz* pfuzz = pit.GetFuzzFrom();
  1321.                     if (pfuzz) {
  1322.                         if (fuzz_from) {
  1323.                             CRef<CInt_fuzz> f(new CInt_fuzz);
  1324.                             f->Assign(*fuzz_from);
  1325.                             f->Subtract(*pfuzz, abs_from, abs_from);
  1326.                             if (f->IsP_m()  &&  !f->GetP_m() ) {
  1327.                                 fuzz_from.Reset(); // cancelled
  1328.                             } else {
  1329.                                 fuzz_from = f;
  1330.                             }
  1331.                         } else {
  1332.                             fuzz_from = pfuzz->Negative(abs_from);
  1333.                         }
  1334.                     }
  1335.                 }
  1336.             } else {
  1337.                 abs_from  = prange.GetFrom();
  1338.                 // fuzz_from = pit.GetFuzzFrom();
  1339.                 CRef<CInt_fuzz> f(new CInt_fuzz);
  1340.                 f->SetLim(CInt_fuzz::eLim_lt);
  1341.                 fuzz_from = f;
  1342.             }
  1343.             if (crange.GetTo() <= prange.GetTo()) {
  1344.                 abs_to  = crange.GetTo();
  1345.                 fuzz_to = cit.GetFuzzTo();
  1346.                 if (abs_to == prange.GetTo()) {
  1347.                     // subtract out parent fuzz, if any
  1348.                     const CInt_fuzz* pfuzz = pit.GetFuzzTo();
  1349.                     if (pfuzz) {
  1350.                         if (fuzz_to) {
  1351.                             CRef<CInt_fuzz> f(new CInt_fuzz);
  1352.                             f->Assign(*fuzz_to);
  1353.                             f->Subtract(*pfuzz, abs_to, abs_to);
  1354.                             if (f->IsP_m()  &&  !f->GetP_m() ) {
  1355.                                 fuzz_to.Reset(); // cancelled
  1356.                             } else {
  1357.                                 fuzz_to = f;
  1358.                             }
  1359.                         } else {
  1360.                             fuzz_to = pfuzz->Negative(abs_to);
  1361.                         }
  1362.                     }
  1363.                 }
  1364.             } else {
  1365.                 abs_to  = prange.GetTo();
  1366.                 // fuzz_to = pit.GetFuzzTo();
  1367.                 CRef<CInt_fuzz> f(new CInt_fuzz);
  1368.                 f->SetLim(CInt_fuzz::eLim_gt);
  1369.                 fuzz_to = f;
  1370.             }
  1371.             if (abs_from <= abs_to) {
  1372.                 if (IsReverse(pstrand)) {
  1373.                     TSeqPos sigma = pos + prange.GetTo();
  1374.                     intersection->SetFrom(sigma - abs_to);
  1375.                     intersection->SetTo  (sigma - abs_from);
  1376.                     if (fuzz_from) {
  1377.                         intersection->SetFuzz_to().AssignTranslated
  1378.                             (*fuzz_from, intersection->GetTo(), abs_from);
  1379.                         intersection->SetFuzz_to().Negate
  1380.                             (intersection->GetTo());
  1381.                     }
  1382.                     if (fuzz_to) {
  1383.                         intersection->SetFuzz_from().AssignTranslated
  1384.                             (*fuzz_to, intersection->GetFrom(), abs_to);
  1385.                         intersection->SetFuzz_from().Negate
  1386.                             (intersection->GetFrom());
  1387.                     }
  1388.                     if (cstrand == eNa_strand_unknown) {
  1389.                         intersection->SetStrand(pstrand);
  1390.                     } else {
  1391.                         intersection->SetStrand(Reverse(cstrand));
  1392.                     }
  1393.                 } else {
  1394.                     TSignedSeqPos delta = pos - prange.GetFrom();
  1395.                     intersection->SetFrom(abs_from + delta);
  1396.                     intersection->SetTo  (abs_to   + delta);
  1397.                     if (fuzz_from) {
  1398.                         intersection->SetFuzz_from().AssignTranslated
  1399.                             (*fuzz_from, intersection->GetFrom(), abs_from);
  1400.                     }
  1401.                     if (fuzz_to) {
  1402.                         intersection->SetFuzz_to().AssignTranslated
  1403.                             (*fuzz_to, intersection->GetTo(), abs_to);
  1404.                     }
  1405.                     if (cstrand == eNa_strand_unknown) {
  1406.                         intersection->SetStrand(pstrand);
  1407.                     } else {
  1408.                         intersection->SetStrand(cstrand);
  1409.                     }
  1410.                 }
  1411.                 // add to m_Ranges, combining with the previous
  1412.                 // interval if possible
  1413.                 if ( !(flags & fNoMerge)  &&  !m_Ranges.empty()
  1414.                     &&  SameOrientation(intersection->GetStrand(),
  1415.                                         m_Ranges.back()->GetStrand()) ) {
  1416.                     if (m_Ranges.back()->GetTo() == intersection->GetFrom() - 1
  1417.                         &&  !IsReverse(intersection->GetStrand()) ) {
  1418.                         m_Ranges.back()->SetTo(intersection->GetTo());
  1419.                         if (intersection->IsSetFuzz_to()) {
  1420.                             m_Ranges.back()->SetFuzz_to
  1421.                                 (intersection->SetFuzz_to());
  1422.                         } else {
  1423.                             m_Ranges.back()->ResetFuzz_to();
  1424.                         }
  1425.                     } else if (m_Ranges.back()->GetFrom()
  1426.                                == intersection->GetTo() + 1
  1427.                                &&  IsReverse(intersection->GetStrand())) {
  1428.                         m_Ranges.back()->SetFrom(intersection->GetFrom());
  1429.                         if (intersection->IsSetFuzz_from()) {
  1430.                             m_Ranges.back()->SetFuzz_from
  1431.                                 (intersection->SetFuzz_from());
  1432.                         } else {
  1433.                             m_Ranges.back()->ResetFuzz_from();
  1434.                         }
  1435.                     } else {
  1436.                         m_Ranges.push_back(intersection);
  1437.                     }
  1438.                 } else {
  1439.                     m_Ranges.push_back(intersection);
  1440.                 }
  1441.             }
  1442.             pos += prange.GetLength();
  1443.         }
  1444.     }
  1445. }
  1448. // Bother trying to merge?
  1449. CRef<CSeq_loc> SRelLoc::Resolve(const CSeq_loc& new_parent, CScope* scope,
  1450.                                 SRelLoc::TFlags /* flags */)
  1451.     const
  1452. {
  1453.     typedef CSeq_loc::TRange TRange0;
  1454.     CRef<CSeq_loc> result(new CSeq_loc);
  1455.     CSeq_loc_mix&  mix = result->SetMix();
  1456.     ITERATE (TRanges, it, m_Ranges) {
  1457.         _ASSERT((*it)->GetFrom() <= (*it)->GetTo());
  1458.         TSeqPos pos = 0, start = (*it)->GetFrom();
  1459.         bool    keep_going = true;
  1460.         for (CSeq_loc_CI pit(new_parent);  pit;  ++pit) {
  1461.             TRange0 prange = pit.GetRange();
  1462.             if (prange.IsWholeTo()  &&  scope) {
  1463.                 // determine actual end
  1464.                 prange.SetToOpen(sequence::GetLength(pit.GetSeq_id(), scope));
  1465.             }
  1466.             TSeqPos length = prange.GetLength();
  1467.             if (start >= pos  &&  start < pos + length) {
  1468.                 TSeqPos              from, to;
  1469.                 CConstRef<CInt_fuzz> fuzz_from, fuzz_to;
  1470.                 ENa_strand           strand;
  1471.                 if (IsReverse(pit.GetStrand())) {
  1472.                     TSeqPos sigma = pos + prange.GetTo();
  1473.                     from = sigma - (*it)->GetTo();
  1474.                     to   = sigma - start;
  1475.                     if (from < prange.GetFrom()  ||  from > sigma) {
  1476.                         from = prange.GetFrom();
  1477.                         keep_going = true;
  1478.                     } else {
  1479.                         keep_going = false;
  1480.                     }
  1481.                     if ( !(*it)->IsSetStrand()
  1482.                         ||  (*it)->GetStrand() == eNa_strand_unknown) {
  1483.                         strand = pit.GetStrand();
  1484.                     } else {
  1485.                         strand = Reverse((*it)->GetStrand());
  1486.                     }
  1487.                     if (from == prange.GetFrom()) {
  1488.                         fuzz_from = pit.GetFuzzFrom();
  1489.                     }
  1490.                     if ( !keep_going  &&  (*it)->IsSetFuzz_to() ) {
  1491.                         CRef<CInt_fuzz> f(new CInt_fuzz);
  1492.                         if (fuzz_from) {
  1493.                             f->Assign(*fuzz_from);
  1494.                         } else {
  1495.                             f->SetP_m(0);
  1496.                         }
  1497.                         f->Subtract((*it)->GetFuzz_to(), from, (*it)->GetTo(),
  1498.                                     CInt_fuzz::eAmplify);
  1499.                         if (f->IsP_m()  &&  !f->GetP_m() ) {
  1500.                             fuzz_from.Reset(); // cancelled
  1501.                         } else {
  1502.                             fuzz_from = f;
  1503.                         }
  1504.                     }
  1505.                     if (to == prange.GetTo()) {
  1506.                         fuzz_to = pit.GetFuzzTo();
  1507.                     }
  1508.                     if (start == (*it)->GetFrom()
  1509.                         &&  (*it)->IsSetFuzz_from()) {
  1510.                         CRef<CInt_fuzz> f(new CInt_fuzz);
  1511.                         if (fuzz_to) {
  1512.                             f->Assign(*fuzz_to);
  1513.                         } else {
  1514.                             f->SetP_m(0);
  1515.                         }
  1516.                         f->Subtract((*it)->GetFuzz_from(), to,
  1517.                                     (*it)->GetFrom(), CInt_fuzz::eAmplify);
  1518.                         if (f->IsP_m()  &&  !f->GetP_m() ) {
  1519.                             fuzz_to.Reset(); // cancelled
  1520.                         } else {
  1521.                             fuzz_to = f;
  1522.                         }
  1523.                     }
  1524.                 } else {
  1525.                     TSignedSeqPos delta = prange.GetFrom() - pos;
  1526.                     from = start          + delta;
  1527.                     to   = (*it)->GetTo() + delta;
  1528.                     if (to > prange.GetTo()) {
  1529.                         to = prange.GetTo();
  1530.                         keep_going = true;
  1531.                     } else {
  1532.                         keep_going = false;
  1533.                     }
  1534.                     if ( !(*it)->IsSetStrand()
  1535.                         ||  (*it)->GetStrand() == eNa_strand_unknown) {
  1536.                         strand = pit.GetStrand();
  1537.                     } else {
  1538.                         strand = (*it)->GetStrand();
  1539.                     }
  1540.                     if (from == prange.GetFrom()) {
  1541.                         fuzz_from = pit.GetFuzzFrom();
  1542.                     }
  1543.                     if (start == (*it)->GetFrom()
  1544.                         &&  (*it)->IsSetFuzz_from()) {
  1545.                         CRef<CInt_fuzz> f(new CInt_fuzz);
  1546.                         if (fuzz_from) {
  1547.                             f->Assign(*fuzz_from);
  1548.                             f->Add((*it)->GetFuzz_from(), from,
  1549.                                    (*it)->GetFrom());
  1550.                         } else {
  1551.                             f->AssignTranslated((*it)->GetFuzz_from(), from,
  1552.                                                 (*it)->GetFrom());
  1553.                         }
  1554.                         if (f->IsP_m()  &&  !f->GetP_m() ) {
  1555.                             fuzz_from.Reset(); // cancelled
  1556.                         } else {
  1557.                             fuzz_from = f;
  1558.                         }
  1559.                     }
  1560.                     if (to == prange.GetTo()) {
  1561.                         fuzz_to = pit.GetFuzzTo();
  1562.                     }
  1563.                     if ( !keep_going  &&  (*it)->IsSetFuzz_to() ) {
  1564.                         CRef<CInt_fuzz> f(new CInt_fuzz);
  1565.                         if (fuzz_to) {
  1566.                             f->Assign(*fuzz_to);
  1567.                             f->Add((*it)->GetFuzz_to(), to, (*it)->GetTo());
  1568.                         } else {
  1569.                             f->AssignTranslated((*it)->GetFuzz_to(), to,
  1570.                                                 (*it)->GetTo());
  1571.                         }
  1572.                         if (f->IsP_m()  &&  !f->GetP_m() ) {
  1573.                             fuzz_to.Reset(); // cancelled
  1574.                         } else {
  1575.                             fuzz_to = f;
  1576.                         }
  1577.                     }
  1578.                 }
  1579.                 if (from == to
  1580.                     &&  (fuzz_from == fuzz_to
  1581.                          ||  (fuzz_from.GetPointer()  &&  fuzz_to.GetPointer()
  1582.                               &&  fuzz_from->Equals(*fuzz_to)))) {
  1583.                     // just a point
  1584.                     CRef<CSeq_loc> loc(new CSeq_loc);
  1585.                     CSeq_point& point = loc->SetPnt();
  1586.                     point.SetPoint(from);
  1587.                     if (strand != eNa_strand_unknown) {
  1588.                         point.SetStrand(strand);
  1589.                     }
  1590.                     if (fuzz_from) {
  1591.                         point.SetFuzz().Assign(*fuzz_from);
  1592.                     }
  1593.                     point.SetId().Assign(pit.GetSeq_id());
  1594.                     mix.Set().push_back(loc);
  1595.                 } else {
  1596.                     CRef<CSeq_loc> loc(new CSeq_loc);
  1597.                     CSeq_interval& ival = loc->SetInt();
  1598.                     ival.SetFrom(from);
  1599.                     ival.SetTo(to);
  1600.                     if (strand != eNa_strand_unknown) {
  1601.                         ival.SetStrand(strand);
  1602.                     }
  1603.                     if (fuzz_from) {
  1604.                         ival.SetFuzz_from().Assign(*fuzz_from);
  1605.                     }
  1606.                     if (fuzz_to) {
  1607.                         ival.SetFuzz_to().Assign(*fuzz_to);
  1608.                     }
  1609.                     ival.SetId().Assign(pit.GetSeq_id());
  1610.                     mix.Set().push_back(loc);
  1611.                 }
  1612.                 if (keep_going) {
  1613.                     start = pos + length;
  1614.                 } else {
  1615.                     break;
  1616.                 }
  1617.             }
  1618.             pos += length;
  1619.         }
  1620.         if (keep_going) {
  1621.             TSeqPos total_length;
  1622.             string  label;
  1623.             new_parent.GetLabel(&label);
  1624.             try {
  1625.                 total_length = sequence::GetLength(new_parent, scope);
  1626.                 ERR_POST(Warning << "SRelLoc::Resolve: Relative position "
  1627.                          << start << " exceeds length (" << total_length
  1628.                          << ") of parent location " << label);
  1629.             } catch (sequence::CNoLength) {
  1630.                 ERR_POST(Warning << "SRelLoc::Resolve: Relative position "
  1631.                          << start
  1632.                          << " exceeds length (???) of parent location "
  1633.                          << label);
  1634.             }            
  1635.         }
  1636.     }
  1637.     // clean up output
  1638.     switch (mix.Get().size()) {
  1639.     case 0:
  1640.         result->SetNull();
  1641.         break;
  1642.     case 1:
  1643.         result.Reset(mix.Set().front());
  1644.         break;
  1645.     default:
  1646.         break;
  1647.     }
  1648.     return result;
  1649. }
  1652. //============================================================================//
  1653. //                                 SeqSearch                                  //
  1654. //============================================================================//
  1656. // Public:
  1657. // =======
  1659. // Constructors and Destructors:
  1660. CSeqSearch::CSeqSearch(IClient *client, bool allow_mismatch) :
  1661.     m_AllowOneMismatch(allow_mismatch),
  1662.     m_MaxPatLen(0),
  1663.     m_Client(client)
  1664. {
  1665.     InitializeMaps();
  1666. }
  1669. CSeqSearch::~CSeqSearch(void) {}
  1672. // Add nucleotide pattern or restriction site to sequence search.
  1673. // Uses ambiguity codes, e.g., R = A and G, H = A, C and T
  1674. void CSeqSearch::AddNucleotidePattern
  1675. (const string& name,
  1676. const string& pat, 
  1677. int cut_site,
  1678. int overhang)
  1679. {
  1680.     string pattern = pat;
  1681.     NStr::TruncateSpaces(pattern);
  1682.     NStr::ToUpper(pattern);
  1684.     // reverse complement pattern to see if it is symetrical
  1685.     string rcomp = ReverseComplement(pattern);
  1687.     bool symmetric = (pattern == rcomp);
  1689.     ENa_strand strand = symmetric ? eNa_strand_both : eNa_strand_plus;
  1691.     AddNucleotidePattern(name, pat, cut_site, overhang, strand);
  1692.     if ( !symmetric ) {
  1693.         AddNucleotidePattern(name, rcomp, pat.length() - cut_site, 
  1694.                              overhang, eNa_strand_minus);
  1695.     }
  1696. }
  1699. // Program passes each character in turn to finite state machine.
  1700. int CSeqSearch::Search
  1701. (int current_state,
  1702.  char ch,
  1703.  int position,
  1704.  int length)
  1705. {
  1706.     if ( !m_Client ) return 0;
  1708.     // on first character, populate state transition table
  1709.     if ( !m_Fsa.IsPrimed() ) {
  1710.         m_Fsa.Prime();
  1711.     }
  1712.     
  1713.     int next_state = m_Fsa.GetNextState(current_state, ch);
  1714.     
  1715.     // report any matches at current state to the client object
  1716.     if ( m_Fsa.IsMatchFound(next_state) ) {
  1717.         ITERATE( vector<CMatchInfo>, it, m_Fsa.GetMatches(next_state) ) {
  1718.   //const CMatchInfo& match = *it;
  1719.             int start = position - it->GetPattern().length() + 1;    
  1721.             // prevent multiple reports of patterns for circular sequences.
  1722.             if ( start < length ) {
  1723.                 m_Client->MatchFound(*it, start);
  1724.             }
  1725.         }
  1726.     }
  1728.     return next_state;
  1729. }
  1732. // Search an entire bioseq.
  1733. void CSeqSearch::Search(const CBioseq_Handle& bsh)
  1734. {
  1735.     if ( !bsh ) return;
  1736.     if ( !m_Client ) return;  // no one to report to, so why search at all.
  1738.     CSeqVector seq_vec = bsh.GetSeqVector(CBioseq_Handle::eCoding_Iupac);
  1739.     size_t seq_len = seq_vec.size();
  1740.     size_t search_len = seq_len;
  1742.     CSeq_inst::ETopology topology =
  1743.         bsh.GetBioseqCore()->GetInst().GetTopology();    
  1744.     if ( topology == CSeq_inst::eTopology_circular ) {
  1745.         search_len += m_MaxPatLen - 1;
  1746.     }
  1747.     
  1748.     int state = m_Fsa.GetInitialState();
  1750.     for ( size_t i = 0; i < search_len; ++i ) {
  1751.         state = Search(state, seq_vec[i % seq_len], i, seq_len);
  1752.     }
  1753. }
  1756. // Private:
  1757. // ========
  1759. // translation finite state machine base codes - ncbi4na
  1760. enum EBaseCode {
  1761.     eBase_gap = 0,
  1762.         eBase_A,      /* A    */
  1763.         eBase_C,      /* C    */
  1764.         eBase_M,      /* AC   */
  1765.         eBase_G,      /* G    */
  1766.         eBase_R,      /* AG   */
  1767.         eBase_S,      /* CG   */
  1768.         eBase_V,      /* ACG  */
  1769.         eBase_T,      /* T    */
  1770.         eBase_W,      /* AT   */
  1771.         eBase_Y,      /* CT   */
  1772.         eBase_H,      /* ACT  */
  1773.         eBase_K,      /* GT   */
  1774.         eBase_D,      /* AGT  */
  1775.         eBase_B,      /* CGT  */
  1776.         eBase_N       /* ACGT */
  1777. };
  1780. map<unsigned char, int> CSeqSearch::sm_CharToEnum;
  1781. map<int, unsigned char> CSeqSearch::sm_EnumToChar;
  1782. map<char, char>         CSeqSearch::sm_Complement;
  1785. void CSeqSearch::InitializeMaps(void)
  1786. {
  1787.     int c, i;
  1788.     static const string iupacna_alphabet      = "-ACMGRSVTWYHKDBN";
  1789.     static const string comp_iupacna_alphabet = "-TGKCYSBAWRDMHVN";
  1791.     if ( sm_CharToEnum.empty() ) {
  1792.         // illegal characters map to eBase_gap (0)
  1793.         for ( c = 0; c < 256; ++c ) {
  1794.             sm_CharToEnum[c] = eBase_gap;
  1795.         }
  1796.         
  1797.         // map iupacna alphabet to EBaseCode
  1798.         for (i = eBase_gap; i <= eBase_N; ++i) {
  1799.             c = iupacna_alphabet[i];
  1800.             sm_CharToEnum[c] = (EBaseCode)i;
  1801.             c = tolower(c);
  1802.             sm_CharToEnum[c] = (EBaseCode)i;
  1803.         }
  1804.         sm_CharToEnum ['U'] = eBase_T;
  1805.         sm_CharToEnum ['u'] = eBase_T;
  1806.         sm_CharToEnum ['X'] = eBase_N;
  1807.         sm_CharToEnum ['x'] = eBase_N;
  1808.         
  1809.         // also map ncbi4na alphabet to EBaseCode
  1810.         for (c = eBase_gap; c <= eBase_N; ++c) {
  1811.             sm_CharToEnum[c] = (EBaseCode)c;
  1812.         }
  1813.     }
  1814.     
  1815.     // map EBaseCode to iupacna alphabet
  1816.     if ( sm_EnumToChar.empty() ) { 
  1817.         for (i = eBase_gap; i <= eBase_N; ++i) {
  1818.             sm_EnumToChar[i] = iupacna_alphabet[i];
  1819.         }
  1820.     }
  1822.     // initialize table to convert character to complement character
  1823.     if ( sm_Complement.empty() ) {
  1824.         int len = iupacna_alphabet.length();
  1825.         for ( i = 0; i < len; ++i ) {
  1826.             sm_Complement.insert(make_pair(iupacna_alphabet[i], comp_iupacna_alphabet[i]));
  1827.         }
  1828.     }
  1829. }
  1832. string CSeqSearch::ReverseComplement(const string& pattern) const
  1833. {
  1834.     size_t len = pattern.length();
  1835.     string rcomp = pattern;
  1836.     rcomp.resize(len);
  1838.     // calculate the complement
  1839.     for ( size_t i = 0; i < len; ++i ) {
  1840.         rcomp[i] = sm_Complement[pattern[i]];
  1841.     }
  1843.     // reverse the complement
  1844.     reverse(rcomp.begin(), rcomp.end());
  1846.     return rcomp;
  1847. }
  1850. void CSeqSearch::AddNucleotidePattern
  1851. (const string& name,
  1852. const string& pattern, 
  1853. int cut_site,
  1854. int overhang,
  1855. ENa_strand strand)
  1856. {
  1857.     size_t pat_len = pattern.length();
  1858.     string temp;
  1859.     temp.resize(pat_len);
  1860.     CMatchInfo info(name,
  1861.                     CNcbiEmptyString::Get(), 
  1862.                     cut_site, 
  1863.                     overhang, 
  1864.                     strand);
  1866.     if ( pat_len > m_MaxPatLen ) {
  1867.         m_MaxPatLen = pat_len;
  1868.     }
  1870.     ExpandPattern(pattern, temp, 0, pat_len, info);
  1871. }
  1874. void CSeqSearch::ExpandPattern
  1875. (const string& pattern,
  1876. string& temp,
  1877. int position,
  1878. int pat_len,
  1879. CMatchInfo& info)
  1880. {
  1881.     static EBaseCode expension[] = { eBase_A, eBase_C, eBase_G, eBase_T };
  1883.     if ( position < pat_len ) {
  1884.         int code = sm_CharToEnum[(int)pattern[position]];
  1885.         
  1886.         for ( int i = 0; i < 4; ++i ) {
  1887.             if ( code & expension[i] ) {
  1888.                 temp[position] = sm_EnumToChar[expension[i]];
  1889.                 ExpandPattern(pattern, temp, position + 1, pat_len, info);
  1890.             }
  1891.         }
  1892.     } else { // recursion base
  1893.         // when position reaches pattern length, store one expanded string.
  1894.         info.m_Pattern = temp;
  1895.         m_Fsa.AddWord(info.m_Pattern, info);
  1897.         if ( m_AllowOneMismatch ) {
  1898.             char ch;
  1899.             // put 'N' at every position if a single mismatch is allowed.
  1900.             for ( int i = 0; i < pat_len; ++i ) {
  1901.                 ch = temp[i];
  1902.                 temp[i] = 'N';
  1904.                 info.m_Pattern = temp;
  1905.                 m_Fsa.AddWord(pattern, info);
  1907.                 // restore proper character, go on to put N in next position.
  1908.                 temp[i] = ch;
  1909.             }
  1910.         }
  1911.     }
  1912. }
  1917. END_SCOPE(objects)
  1918. END_NCBI_SCOPE
  1920. /*
  1921. * ===========================================================================
  1922. * $Log: sequence.cpp,v $
  1923. * Revision 1000.3  2004/06/01 19:25:30  gouriano
  1924. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.82
  1925. *
  1926. * Revision 1.82  2004/05/27 13:48:24  jcherry
  1927. * Replaced some calls to deprecated CBioseq_Handle::GetBioseq() with
  1928. * GetCompleteBioseq() or GetBioseqCore()
  1929. *
  1930. * Revision 1.81  2004/05/25 21:06:34  johnson
  1931. * Bug fix in x_Translate
  1932. *
  1933. * Revision 1.80  2004/05/21 21:42:14  gorelenk
  1934. * Added PCH ncbi_pch.hpp
  1935. *
  1936. * Revision 1.79  2004/05/06 17:39:43  shomrat
  1937. * Fixes to SeqLocMerge
  1938. *
  1939. * Revision 1.78  2004/04/06 14:03:15  dicuccio
  1940. * Added API to extract the single Org-ref from a bioseq handle.  Added API to
  1941. * retrieve a single tax-id from a bioseq handle
  1942. *
  1943. * Revision 1.77  2004/04/05 15:56:14  grichenk
  1944. * Redesigned CAnnotTypes_CI: moved all data and data collecting
  1945. * functions to CAnnotDataCollector. CAnnotTypes_CI is no more
  1946. * inherited from SAnnotSelector.
  1947. *
  1948. * Revision 1.76  2004/03/25 20:02:30  vasilche
  1949. * Added several method variants with CBioseq_Handle as argument.
  1950. *
  1951. * Revision 1.75  2004/03/01 18:28:18  dicuccio
  1952. * Changed sequence::Compare() such that a seq-interval of length 1 and a
  1953. * corresponding seq-point compare as the same
  1954. *
  1955. * Revision 1.74  2004/03/01 18:24:22  shomrat
  1956. * Removed branching in the main cdregion translation loop; Added alternative start flag
  1957. *
  1958. * Revision 1.73  2004/02/19 18:16:48  shomrat
  1959. * Implemented SeqlocRevCmp
  1960. *
  1961. * Revision 1.72  2004/02/09 14:43:18  vasilche
  1962. * Added missing typename keyword.
  1963. *
  1964. * Revision 1.71  2004/02/05 19:41:46  shomrat
  1965. * Convenience functions for popular overlapping types
  1966. *
  1967. * Revision 1.70  2004/02/04 18:05:41  grichenk
  1968. * Added annotation filtering by set of types/subtypes.
  1969. * Renamed *Choice to *Type in SAnnotSelector.
  1970. *
  1971. * Revision 1.69  2004/01/30 17:49:29  ucko
  1972. * Add missing "typename"
  1973. *
  1974. * Revision 1.68  2004/01/30 17:22:53  dicuccio
  1975. * Added sequence::GetId(const CBioseq_Handle&) - returns a selected ID class
  1976. * (best, GI).  Added CSeqTranslator - utility class to translate raw sequence
  1977. * data
  1978. *
  1979. * Revision 1.67  2004/01/28 17:19:04  shomrat
  1980. * Implemented SeqLocMerge
  1981. *
  1982. * Revision 1.66  2003/12/16 19:37:43  shomrat
  1983. * Retrieve encoding feature and bioseq of a protein
  1984. *
  1985. * Revision 1.65  2003/11/19 22:18:05  grichenk
  1986. * All exceptions are now CException-derived. Catch "exception" rather
  1987. * than "runtime_error".
  1988. *
  1989. * Revision 1.64  2003/10/17 20:55:27  ucko
  1990. * SRelLoc::Resolve: fix a fuzz-handling paste-o.
  1991. *
  1992. * Revision 1.63  2003/10/16 11:55:19  dicuccio
  1993. * Fix for brain-dead MSVC and ambiguous operator&&
  1994. *
  1995. * Revision 1.62  2003/10/15 19:52:18  ucko
  1996. * More adjustments to SRelLoc: support fuzz, opposite-strand children,
  1997. * and resolving against an alternate parent.
  1998. *
  1999. * Revision 1.61  2003/10/08 21:08:38  ucko
  2000. * CCdregion_translate: take const Bioseq_Handles, since there's no need
  2001. * to modify them.
  2002. * TestForOverlap: don't consider sequences circular if
  2003. * circular_len == kInvalidSeqPos
  2004. *
  2005. * Revision 1.60  2003/09/22 18:38:14  grichenk
  2006. * Fixed circular seq-locs processing by TestForOverlap()
  2007. *
  2008. * Revision 1.59  2003/07/17 21:00:50  vasilche
  2009. * Added missing include <list>
  2010. *
  2011. * Revision 1.58  2003/06/19 17:11:43  ucko
  2012. * SRelLoc::SRelLoc: remember to update the base position when skipping
  2013. * parent ranges whose IDs don't match.
  2014. *
  2015. * Revision 1.57  2003/06/13 17:23:32  grichenk
  2016. * Added special processing of multi-ID seq-locs in TestForOverlap()
  2017. *
  2018. * Revision 1.56  2003/06/05 18:08:36  shomrat
  2019. * Allow empty location when computing SeqLocPartial; Adjust GetSeqData in cdregion translation
  2020. *
  2021. * Revision 1.55  2003/06/02 18:58:25  dicuccio
  2022. * Fixed #include directives
  2023. *
  2024. * Revision 1.54  2003/06/02 18:53:32  dicuccio
  2025. * Fixed bungled commit
  2026. *
  2027. * Revision 1.52  2003/05/27 19:44:10  grichenk
  2028. * Added CSeqVector_CI class
  2029. *
  2030. * Revision 1.51  2003/05/15 19:27:02  shomrat
  2031. * Compare handle only if both valid; Check IsLim before GetLim
  2032. *
  2033. * Revision 1.50  2003/05/09 15:37:00  ucko
  2034. * Take const CBioseq_Handle references in CFastaOstream::Write et al.
  2035. *
  2036. * Revision 1.49  2003/05/06 19:34:36  grichenk
  2037. * Fixed GetStrand() (worked fine only for plus and unknown strands)
  2038. *
  2039. * Revision 1.48  2003/05/01 17:55:17  ucko
  2040. * Fix GetLength(const CSeq_id&, CScope*) to return ...::max() rather
  2041. * than throwing if it can't resolve the ID to a handle.
  2042. *
  2043. * Revision 1.47  2003/04/24 16:15:58  vasilche
  2044. * Added missing includes and forward class declarations.
  2045. *
  2046. * Revision 1.46  2003/04/16 19:44:26  grichenk
  2047. * More fixes to TestForOverlap() and GetStrand()
  2048. *
  2049. * Revision 1.45  2003/04/15 20:11:21  grichenk
  2050. * Fixed strands problem in TestForOverlap(), replaced type
  2051. * iterator with container iterator in s_GetStrand().
  2052. *
  2053. * Revision 1.44  2003/04/03 19:03:17  grichenk
  2054. * Two more cases to revert locations in GetBestOverlappingFeat()
  2055. *
  2056. * Revision 1.43  2003/04/02 16:58:59  grichenk
  2057. * Change location and feature in GetBestOverlappingFeat()
  2058. * for eOverlap_Subset.
  2059. *
  2060. * Revision 1.42  2003/03/25 22:00:20  grichenk
  2061. * Added strand checking to TestForOverlap()
  2062. *
  2063. * Revision 1.41  2003/03/18 21:48:35  grichenk
  2064. * Removed obsolete class CAnnot_CI
  2065. *
  2066. * Revision 1.40  2003/03/11 16:00:58  kuznets
  2067. * iterate -> ITERATE
  2068. *
  2069. * Revision 1.39  2003/02/19 16:25:14  grichenk
  2070. * Check strands in GetBestOverlappingFeat()
  2071. *
  2072. * Revision 1.38  2003/02/14 15:41:00  shomrat
  2073. * Minor implementation changes in SeqLocPartialTest
  2074. *
  2075. * Revision 1.37  2003/02/13 14:35:40  grichenk
  2076. * + eOverlap_Contains
  2077. *
  2078. * Revision 1.36  2003/02/10 15:54:01  grichenk
  2079. * Use CFeat_CI->GetMappedFeature() and GetOriginalFeature()
  2080. *
  2081. * Revision 1.35  2003/02/06 22:26:27  vasilche
  2082. * Use CSeq_id::Assign().
  2083. *
  2084. * Revision 1.34  2003/02/06 20:59:16  shomrat
  2085. * Bug fix in SeqLocPartialCheck
  2086. *
  2087. * Revision 1.33  2003/01/22 21:02:23  ucko
  2088. * Fix stupid logic error in SRelLoc's constructor; change LocationOffset
  2089. * to return -1 rather than crashing if the locations don't intersect.
  2090. *
  2091. * Revision 1.32  2003/01/22 20:15:02  vasilche
  2092. * Removed compiler warning.
  2093. *
  2094. * Revision 1.31  2003/01/22 18:17:09  ucko
  2095. * SRelLoc::SRelLoc: change intersection to a CRef, so we don't have to
  2096. * worry about it going out of scope while still referenced (by m_Ranges).
  2097. *
  2098. * Revision 1.30  2003/01/08 20:43:10  ucko
  2099. * Adjust SRelLoc to use (ID-less) Seq-intervals for ranges, so that it
  2100. * will be possible to add support for fuzz and strandedness/orientation.
  2101. *
  2102. * Revision 1.29  2002/12/30 19:38:35  vasilche
  2103. * Optimized CGenbankWriter::WriteSequence.
  2104. * Implemented GetBestOverlappingFeat() with CSeqFeatData::ESubtype selector.
  2105. *
  2106. * Revision 1.28  2002/12/26 21:45:29  grichenk
  2107. * + GetBestOverlappingFeat()
  2108. *
  2109. * Revision 1.27  2002/12/26 21:17:06  dicuccio
  2110. * Minor tweaks to avoid compiler warnings in MSVC (remove unused variables)
  2111. *
  2112. * Revision 1.26  2002/12/20 17:14:18  kans
  2113. * added SeqLocPartialCheck
  2114. *
  2115. * Revision 1.25  2002/12/19 20:24:55  grichenk
  2116. * Updated usage of CRange<>
  2117. *
  2118. * Revision 1.24  2002/12/10 16:56:41  ucko
  2119. * CFastaOstream::WriteTitle: restore leading > accidentally dropped in R1.19.
  2120. *
  2121. * Revision 1.23  2002/12/10 16:34:45  kans
  2122. * implement code break processing in TranslateCdregion
  2123. *
  2124. * Revision 1.22  2002/12/09 20:38:41  ucko
  2125. * +sequence::LocationOffset
  2126. *
  2127. * Revision 1.21  2002/12/06 15:36:05  grichenk
  2128. * Added overlap type for annot-iterators
  2129. *
  2130. * Revision 1.20  2002/12/04 15:38:22  ucko
  2131. * SourceToProduct, ProductToSource: just check whether the feature is a coding
  2132. * region rather than attempting to determine molecule types; drop s_DeduceMol.
  2133. *
  2134. * Revision 1.19  2002/11/26 15:14:04  dicuccio
  2135. * Changed CFastaOStream::WriteTitle() to make use of CSeq_id::GetStringDescr().
  2136. *
  2137. * Revision 1.18  2002/11/25 21:24:46  grichenk
  2138. * Added TestForOverlap() function.
  2139. *
  2140. * Revision 1.17  2002/11/18 19:59:23  shomrat
  2141. * Add CSeqSearch - a nucleotide search utility
  2142. *
  2143. * Revision 1.16  2002/11/12 20:00:25  ucko
  2144. * +SourceToProduct, ProductToSource, SRelLoc
  2145. *
  2146. * Revision 1.15  2002/10/23 19:22:39  ucko
  2147. * Move the FASTA reader from objects/util/sequence.?pp to
  2148. * objects/seqset/Seq_entry.?pp because it doesn't need the OM.
  2149. *
  2150. * Revision 1.14  2002/10/23 18:23:59  ucko
  2151. * Clean up #includes.
  2152. * Add a FASTA reader (known to compile, but not otherwise tested -- take care)
  2153. *
  2154. * Revision 1.13  2002/10/23 16:41:12  clausen
  2155. * Fixed warning in WorkShop53
  2156. *
  2157. * Revision 1.12  2002/10/23 15:33:50  clausen
  2158. * Fixed local variable scope warnings
  2159. *
  2160. * Revision 1.11  2002/10/08 12:35:37  clausen
  2161. * Fixed bugs in GetStrand(), ChangeSeqId() & SeqLocCheck()
  2162. *
  2163. * Revision 1.10  2002/10/07 17:11:16  ucko
  2164. * Fix usage of ERR_POST (caught by KCC)
  2165. *
  2166. * Revision 1.9  2002/10/03 18:44:09  clausen
  2167. * Removed extra whitespace
  2168. *
  2169. * Revision 1.8  2002/10/03 16:33:55  clausen
  2170. * Added functions needed by validate
  2171. *
  2172. * Revision 1.7  2002/09/13 15:30:43  ucko
  2173. * Change resize(0) to erase() at Denis's suggestion.
  2174. *
  2175. * Revision 1.6  2002/09/13 14:28:34  ucko
  2176. * string::clear() doesn't exist under g++ 2.95, so use resize(0) instead. :-/
  2177. *
  2178. * Revision 1.5  2002/09/12 21:39:13  kans
  2179. * added CCdregion_translate and CCdregion_translate
  2180. *
  2181. * Revision 1.4  2002/09/03 21:27:04  grichenk
  2182. * Replaced bool arguments in CSeqVector constructor and getters
  2183. * with enums.
  2184. *
  2185. * Revision 1.3  2002/08/27 21:41:15  ucko
  2186. * Add CFastaOstream.
  2187. *
  2188. * Revision 1.2  2002/06/07 16:11:09  ucko
  2189. * Move everything into the "sequence" namespace.
  2190. * Drop the anonymous-namespace business; "sequence" should provide
  2191. * enough protection, and anonymous namespaces ironically interact poorly
  2192. * with WorkShop's caching when rebuilding shared libraries.
  2193. *
  2194. * Revision 1.1  2002/06/06 18:41:41  clausen
  2195. * Initial version
  2196. *
  2197. * ===========================================================================
  2198. */