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

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

生物技术

开发平台：

C/C++

seq_map.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: seq_map.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 19:24:15  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.53
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: seq_map.cpp,v 1000.4 2004/06/01 19:24:15 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors: Aleksey Grichenko, Michael Kimelman, Eugene Vasilchenko,
  36. *          Andrei Gourianov
  37. *
  38. * File Description:
  39. *   Sequence map for the Object Manager. Describes sequence as a set of
  40. *   segments of different types (data, reference, gap or end).
  41. *
  42. */
  44. #include <ncbi_pch.hpp>
  45. #include <objmgr/seq_map.hpp>
  46. #include <objmgr/seq_map_ci.hpp>
  47. #include <objmgr/seq_map_ext.hpp>
  48. #include <objmgr/seq_id_mapper.hpp>
  49. #include <objmgr/scope.hpp>
  50. #include <objmgr/bioseq_handle.hpp>
  52. #include <objects/seq/Bioseq.hpp>
  53. #include <objects/seq/Seq_data.hpp>
  54. #include <objects/seq/Seq_inst.hpp>
  55. #include <objects/seq/Delta_ext.hpp>
  56. #include <objects/seq/Delta_seq.hpp>
  57. #include <objects/seq/Seq_literal.hpp>
  58. #include <objects/seq/Seq_ext.hpp>
  59. #include <objects/seq/Seg_ext.hpp>
  60. #include <objects/seq/Ref_ext.hpp>
  62. #include <objects/seqloc/Seq_loc.hpp>
  63. #include <objects/seqloc/Seq_point.hpp>
  64. #include <objects/seqloc/Seq_loc_mix.hpp>
  65. #include <objects/seqloc/Seq_loc_equiv.hpp>
  66. #include <objects/seqloc/Seq_interval.hpp>
  67. #include <objects/seqloc/Packed_seqint.hpp>
  68. #include <objects/seqloc/Packed_seqpnt.hpp>
  70. #include <algorithm>
  72. BEGIN_NCBI_SCOPE
  73. BEGIN_SCOPE(objects)
  75. ////////////////////////////////////////////////////////////////////
  76. //  CSeqMap::CSegment
  78. inline
  79. CSeqMap::CSegment::CSegment(ESegmentType seg_type,
  80.                             TSeqPos length)
  81.     : m_Position(kInvalidSeqPos),
  82.       m_Length(length),
  83.       m_SegType(seg_type),
  84.       m_RefMinusStrand(false),
  85.       m_RefPosition(0)
  86. {
  87. }
  89. ////////////////////////////////////////////////////////////////////
  90. //  CSeqMap
  93. CSeqMap::CSeqMap(void)
  94.     : m_Resolved(0),
  95.       m_Mol(CSeq_inst::eMol_not_set),
  96.       m_SeqLength(kInvalidSeqPos)
  97. {
  98. }
  101. CSeqMap::CSeqMap(CSeqMap* /*parent*/, size_t /*index*/)
  102.     : m_Resolved(0),
  103.       m_Mol(CSeq_inst::eMol_not_set),
  104.       m_SeqLength(kInvalidSeqPos)
  105. {
  106. }
  109. CSeqMap::CSeqMap(const CSeq_loc& ref)
  110.     : m_Resolved(0),
  111.       m_Mol(CSeq_inst::eMol_not_set),
  112.       m_SeqLength(kInvalidSeqPos)
  113. {
  114.     x_AddEnd();
  115.     x_Add(ref);
  116.     x_AddEnd();
  117. }
  120. CSeqMap::CSeqMap(const CSeq_data& data, TSeqPos length)
  121.     : m_Resolved(0),
  122.       m_Mol(CSeq_inst::eMol_not_set),
  123.       m_SeqLength(kInvalidSeqPos)
  124. {
  125.     x_AddEnd();
  126.     x_Add(data, length);
  127.     x_AddEnd();
  128. }
  131. CSeqMap::CSeqMap(TSeqPos length)
  132.     : m_Resolved(0),
  133.       m_Mol(CSeq_inst::eMol_not_set),
  134.       m_SeqLength(length)
  135. {
  136.     x_AddEnd();
  137.     x_AddGap(length);
  138.     x_AddEnd();
  139. }
  142. CSeqMap::~CSeqMap(void)
  143. {
  144.     m_Resolved = 0;
  145.     m_Segments.clear();
  146. }
  149. void CSeqMap::x_GetSegmentException(size_t /*index*/) const
  150. {
  151.     NCBI_THROW(CSeqMapException, eInvalidIndex,
  152.                "Invalid segment index");
  153. }
  156. CSeqMap::CSegment& CSeqMap::x_SetSegment(size_t index)
  157. {
  158.     if ( index >= x_GetSegmentsCount() ) {
  159.         NCBI_THROW(CSeqMapException, eInvalidIndex,
  160.                    "Invalid segment index");
  161.     }
  162.     return m_Segments[index];
  163. }
  166. CBioseq_Handle CSeqMap::x_GetBioseqHandle(const CSegment& seg, CScope* scope) const
  167. {
  168.     if ( !scope ) {
  169.         NCBI_THROW(CSeqMapException, eNullPointer,
  170.                    "Null scope pointer");
  171.     }
  172.     CBioseq_Handle bh = scope->GetBioseqHandle(x_GetRefSeqid(seg));
  173.     if ( !bh ) {
  174.         NCBI_THROW(CSeqMapException, eFail,
  175.                    "Cannot resolve reference");
  176.     }
  177.     return bh;
  178. }
  181. TSeqPos CSeqMap::x_ResolveSegmentLength(size_t index, CScope* scope) const
  182. {
  183.     const CSegment& seg = x_GetSegment(index);
  184.     TSeqPos length = seg.m_Length;
  185.     if ( length == kInvalidSeqPos ) {
  186.         if ( seg.m_SegType == eSeqSubMap ) {
  187.             length = x_GetSubSeqMap(seg, scope)->GetLength(scope);
  188.         }
  189.         else if ( seg.m_SegType == eSeqRef ) {
  190.             length = x_GetBioseqHandle(seg, scope).GetBioseqLength();
  191.         }
  192.         _ASSERT(length != kInvalidSeqPos);
  193.         seg.m_Length = length;
  194.     }
  195.     return length;
  196. }
  199. TSeqPos CSeqMap::x_ResolveSegmentPosition(size_t index, CScope* scope) const
  200. {
  201.     if ( index > x_GetSegmentsCount() ) {
  202.         NCBI_THROW(CSeqMapException, eInvalidIndex,
  203.                    "Invalid segment index");
  204.     }
  205.     size_t resolved = m_Resolved;
  206.     if ( index <= resolved )
  207.         return x_GetSegment(index).m_Position;
  208.     TSeqPos resolved_pos = x_GetSegment(resolved).m_Position;
  209.     do {
  210.         resolved_pos += x_GetSegmentLength(resolved, scope);
  211.         m_Segments[++resolved].m_Position = resolved_pos;
  212.     } while ( resolved < index );
  213.     {{
  214.         CFastMutexGuard guard(m_SeqMap_Mtx);
  215.         if ( m_Resolved < resolved )
  216.             m_Resolved = resolved;
  217.     }}
  218.     return resolved_pos;
  219. }
  222. size_t CSeqMap::x_FindSegment(TSeqPos pos, CScope* scope) const
  223. {
  224.     size_t resolved = m_Resolved;
  225.     TSeqPos resolved_pos = x_GetSegment(resolved).m_Position;
  226.     if ( resolved_pos <= pos ) {
  227.         do {
  228.             if ( resolved >= x_GetSegmentsCount() ) {
  229.                 // end of segments
  230.                 m_Resolved = resolved;
  231.                 return size_t(-1);
  232.             }
  233.             resolved_pos += x_GetSegmentLength(resolved, scope);
  234.             m_Segments[++resolved].m_Position = resolved_pos;
  235.         } while ( resolved_pos <= pos );
  236.         {{
  237.             CFastMutexGuard guard(m_SeqMap_Mtx);
  238.             if ( m_Resolved < resolved )
  239.                 m_Resolved = resolved;
  240.         }}
  241.         return resolved - 1;
  242.     }
  243.     else {
  244.         TSegments::const_iterator end = m_Segments.begin()+resolved;
  245.         TSegments::const_iterator it = 
  246.             upper_bound(m_Segments.begin(), end,
  247.                         pos, SPosLessSegment());
  248.         if ( it == end ) {
  249.             return size_t(-1);
  250.         }
  251.         return it - m_Segments.begin();
  252.     }
  253. }
  256. void CSeqMap::x_LoadObject(const CSegment& seg) const
  257. {
  258.     _ASSERT(seg.m_Position != kInvalidSeqPos);
  259.     if ( !seg.m_RefObject ) {
  260.         NCBI_THROW(CSeqMapException, eFail, "Cannot load part of seqmap");
  261.     }
  262. }
  265. CConstRef<CSeqMap> CSeqMap::x_GetSubSeqMap(const CSegment& seg, CScope* scope,
  266.                                            bool resolveExternal) const
  267. {
  268.     CConstRef<CSeqMap> ret;
  269.     if ( seg.m_SegType == eSeqSubMap ) {
  270.         if ( !seg.m_RefObject ) {
  271.             x_LoadObject(seg);
  272.         }
  273.         if ( !seg.m_RefObject ) {
  274.         NCBI_THROW(CSeqMapException, eNullPointer, "Null CSeqMap pointer");
  275.         }
  276.         ret.Reset(static_cast<const CSeqMap*>(seg.m_RefObject.GetPointer()));
  277.     }
  278.     else if ( resolveExternal && seg.m_SegType == eSeqRef ) {
  279.         ret.Reset(&x_GetBioseqHandle(seg, scope).GetSeqMap());
  280.     }
  281.     return ret;
  282. }
  285. void CSeqMap::x_SetSubSeqMap(size_t /*index*/, CSeqMap_Delta_seqs* /*subMap*/)
  286. {
  287.     // not valid in generic seq map -> incompatible objects
  288.     NCBI_THROW(CSeqMapException, eDataError, "Invalid parent map");
  289. }
  292. const CSeq_data& CSeqMap::x_GetSeq_data(const CSegment& seg) const
  293. {
  294.     if ( seg.m_SegType == eSeqData ) {
  295.         if ( !seg.m_RefObject ) {
  296.             x_LoadObject(seg);
  297.         }
  298.         if ( !seg.m_RefObject ) {
  299.             NCBI_THROW(CSeqMapException, eNullPointer,
  300.                        "Null CSeq_data pointer");
  301.         }
  302.         return *static_cast<const CSeq_data*>(seg.m_RefObject.GetPointer());
  303.     }
  304.     NCBI_THROW(CSeqMapException, eSegmentTypeError,
  305.                "Invalid segment type");
  306. }
  309. void CSeqMap::x_SetSeq_data(size_t index, CSeq_data& data)
  310. {
  311.     // check segment type
  312.     CSegment& seg = x_SetSegment(index);
  313.     if ( seg.m_SegType != eSeqData ) {
  314.         NCBI_THROW(CSeqMapException, eSegmentTypeError,
  315.                    "Invalid segment type");
  316.     }
  317.     // lock for object modification
  318.     CFastMutexGuard guard(m_SeqMap_Mtx);
  319.     // check for object
  320.     if ( seg.m_RefObject ) {
  321.         NCBI_THROW(CSeqMapException, eDataError, "CSeq_data already set");
  322.     }
  323.     // set object
  324.     seg.m_RefObject.Reset(&data);
  325. }
  328. const CSeq_id& CSeqMap::x_GetRefSeqid(const CSegment& seg) const
  329. {
  330.     if ( seg.m_SegType == eSeqRef ) {
  331.         if ( !seg.m_RefObject ) {
  332.             NCBI_THROW(CSeqMapException, eNullPointer, "Null CSeq_id pointer");
  333.         }
  334.         return static_cast<const CSeq_id&>(*seg.m_RefObject);
  335.     }
  336.     NCBI_THROW(CSeqMapException, eSegmentTypeError,
  337.                "Invalid segment type");
  338. }
  341. TSeqPos CSeqMap::x_GetRefPosition(const CSegment& seg) const
  342. {
  343.     return seg.m_RefPosition;
  344. }
  347. bool CSeqMap::x_GetRefMinusStrand(const CSegment& seg) const
  348. {
  349.     return seg.m_RefMinusStrand;
  350. }
  353. CSeqMap::TSegment_CI CSeqMap::Begin(CScope* scope) const
  354. {
  355.     return TSegment_CI(CConstRef<CSeqMap>(this), scope, 0);
  356. }
  359. CSeqMap::TSegment_CI CSeqMap::End(CScope* scope) const
  360. {
  361.     return TSegment_CI(CConstRef<CSeqMap>(this), scope, GetLength(scope));
  362. }
  365. CSeqMap::TSegment_CI CSeqMap::FindSegment(TSeqPos pos, CScope* scope) const
  366. {
  367.     return TSegment_CI(CConstRef<CSeqMap>(this), scope, pos);
  368. }
  371. CSeqMap::const_iterator CSeqMap::begin(CScope* scope) const
  372. {
  373.     return Begin(scope);
  374. }
  377. CSeqMap::const_iterator CSeqMap::end(CScope* scope) const
  378. {
  379.     return End(scope);
  380. }
  383. CSeqMap::const_iterator CSeqMap::find(TSeqPos pos, CScope* scope) const
  384. {
  385.     return FindSegment(pos, scope);
  386. }
  389. CSeqMap::TSegment_CI CSeqMap::BeginResolved(CScope* scope, size_t maxResolveCount,
  390.                                             TFlags flags) const
  391. {
  392.     return begin_resolved(scope, maxResolveCount, flags);
  393. }
  396. CSeqMap::TSegment_CI CSeqMap::EndResolved(CScope* scope, size_t maxResolveCount,
  397.                                           TFlags flags) const
  398. {
  399.     return end_resolved(scope, maxResolveCount, flags);
  400. }
  403. CSeqMap::TSegment_CI CSeqMap::FindResolved(TSeqPos pos, CScope* scope,
  404.                                            size_t maxResolveCount,
  405.                                            TFlags flags) const
  406. {
  407.     return find_resolved(pos, scope, maxResolveCount, flags);
  408. }
  411. CSeqMap::TSegment_CI CSeqMap::FindResolved(TSeqPos pos, CScope* scope,
  412.                                            ENa_strand strand,
  413.                                            size_t maxResolveCount,
  414.                                            TFlags flags) const
  415. {
  416.     return find_resolved(pos, scope, strand, maxResolveCount, flags);
  417. }
  420. CSeqMap::const_iterator CSeqMap::begin_resolved(CScope* scope,
  421.                                                 size_t maxResolveCount,
  422.                                                 TFlags flags) const
  423. {
  424.     return TSegment_CI(CConstRef<CSeqMap>(this), scope, 0,
  425.                        maxResolveCount, flags);
  426. }
  429. CSeqMap::const_iterator CSeqMap::end_resolved(CScope* scope,
  430.                                               size_t maxResolveCount,
  431.                                               TFlags flags) const
  432. {
  433.     return TSegment_CI(CConstRef<CSeqMap>(this), scope,
  434.                        GetLength(scope),
  435.                        maxResolveCount, flags);
  436. }
  439. CSeqMap::const_iterator CSeqMap::find_resolved(TSeqPos pos, CScope* scope,
  440.                                                size_t maxResolveCount,
  441.                                                TFlags flags) const
  442. {
  443.     return TSegment_CI(CConstRef<CSeqMap>(this), scope, pos,
  444.                        maxResolveCount, flags);
  445. }
  448. CSeqMap::const_iterator CSeqMap::find_resolved(TSeqPos pos, CScope* scope,
  449.                                                ENa_strand strand,
  450.                                                size_t maxResolveCount,
  451.                                                TFlags flags) const
  452. {
  453.     return TSegment_CI(CConstRef<CSeqMap>(this), scope,
  454.                        pos, strand,
  455.                        maxResolveCount, flags);
  456. }
  459. CSeqMap::TSegment_CI
  460. CSeqMap::ResolvedRangeIterator(CScope* scope,
  461.                                ENa_strand strand,
  462.                                TSeqPos from,
  463.                                TSeqPos length,
  464.                                size_t maxResolveCount,
  465.                                TFlags flags) const
  466. {
  467.     if ( IsReverse(strand) ) {
  468.         from = GetLength(scope) - from - length;
  469.     }
  470.     return TSegment_CI(CConstRef<CSeqMap>(this), scope,
  471.                        SSeqMapSelector()
  472.                        .SetRange(from, length)
  473.                        .SetResolveCount(maxResolveCount)
  474.                        .SetFlags(flags),
  475.                        strand);
  476. }
  479. CSeqMap::TSegment_CI
  480. CSeqMap::ResolvedRangeIterator(CScope* scope,
  481.                                TSeqPos from,
  482.                                TSeqPos length,
  483.                                ENa_strand strand,
  484.                                size_t maxResolveCount,
  485.                                TFlags flags) const
  486. {
  487.     return TSegment_CI(CConstRef<CSeqMap>(this), scope,
  488.                        SSeqMapSelector()
  489.                        .SetRange(from, length)
  490.                        .SetResolveCount(maxResolveCount)
  491.                        .SetFlags(flags),
  492.                        strand);
  493. }
  496. bool CSeqMap::CanResolveRange(CScope* scope,
  497.                               TSeqPos from,
  498.                               TSeqPos length,
  499.                               ENa_strand strand) const
  500. {
  501.     try {
  502.         TSegment_CI seg(CConstRef<CSeqMap>(this), scope,
  503.                         SSeqMapSelector()
  504.                         .SetRange(from, length)
  505.                         .SetResolveCount(size_t(-1))
  506.                         .SetFlags(fDefaultFlags),
  507.                         strand);
  508.         for ( ; seg; ++seg);
  509.     }
  510.     catch (exception) {
  511.         return false;
  512.     }
  513.     return true;
  514. }
  517. void CSeqMap::DebugDump(CDebugDumpContext /*ddc*/,
  518.                         unsigned int /*depth*/) const
  519. {
  520. #if 0
  521.     ddc.SetFrame("CSeqMap");
  522.     CObject::DebugDump( ddc, depth);
  524.     DebugDumpValue(ddc, "m_FirstUnresolvedPos", m_FirstUnresolvedPos);
  525.     if (depth == 0) {
  526.         DebugDumpValue(ddc, "m_Data.size()", m_Data.size());
  527.     } else {
  528.         DebugDumpValue(ddc, "m_Data.type",
  529.             "vector<CRef<CSegment>>");
  530.         DebugDumpRangeCRef(ddc, "m_Data",
  531.             m_Data.begin(), m_Data.end(), depth);
  532.     }
  533. #endif
  534. }
  537. CConstRef<CSeqMap> CSeqMap::CreateSeqMapForBioseq(const CBioseq& seq)
  538. {
  539.     CConstRef<CSeqMap> ret;
  540.     const CSeq_inst& inst = seq.GetInst();
  541.     if ( inst.IsSetSeq_data() ) {
  542.         ret.Reset(new CSeqMap(inst.GetSeq_data(),
  543.                               inst.GetLength()));
  544.     }
  545.     else if ( inst.IsSetExt() ) {
  546.         const CSeq_ext& ext = inst.GetExt();
  547.         switch (ext.Which()) {
  548.         case CSeq_ext::e_Seg:
  549.             ret.Reset(new CSeqMap_Seq_locs(ext.GetSeg(),
  550.                                            ext.GetSeg().Get()));
  551.             break;
  552.         case CSeq_ext::e_Ref:
  553.             ret.Reset(new CSeqMap(ext.GetRef()));
  554.             break;
  555.         case CSeq_ext::e_Delta:
  556.             ret.Reset(new CSeqMap_Delta_seqs(ext.GetDelta()));
  557.             break;
  558.         case CSeq_ext::e_Map:
  559.             //### Not implemented
  560.             NCBI_THROW(CSeqMapException, eUnimplemented,
  561.                        "CSeq_ext::e_Map -- not implemented");
  562.         default:
  563.             //### Not implemented
  564.             NCBI_THROW(CSeqMapException, eUnimplemented,
  565.                        "CSeq_ext::??? -- not implemented");
  566.         }
  567.     }
  568.     else if ( inst.GetRepr() == CSeq_inst::eRepr_virtual ) {
  569.         // Virtual sequence -- no data, no segments
  570.         // The total sequence is gap
  571.         ret.Reset(new CSeqMap(inst.GetLength()));
  572.     }
  573.     else {
  574.         if ( inst.GetRepr() != CSeq_inst::eRepr_not_set ) {
  575.             NCBI_THROW(CSeqMapException, eDataError,
  576.                        "CSeq_inst.repr of sequence without data "
  577.                        "should be not_set");
  578.         }
  579.         if ( inst.IsSetLength() && inst.GetLength() != 0 ) {
  580.             NCBI_THROW(CSeqMapException, eDataError,
  581.                        "CSeq_inst.length of sequence without data "
  582.                        "should be 0");
  583.         }
  584.         ret.Reset(new CSeqMap(TSeqPos(0)));
  585.     }
  586.     const_cast<CSeqMap&>(*ret).m_Mol = inst.GetMol();
  587.     if ( inst.IsSetLength() ) {
  588.         const_cast<CSeqMap&>(*ret).m_SeqLength = inst.GetLength();
  589.     }
  590.     return ret;
  591. }
  594. CConstRef<CSeqMap> CSeqMap::CreateSeqMapForSeq_loc(const CSeq_loc& loc,
  595.                                                    CScope* scope)
  596. {
  597.     CConstRef<CSeqMap> ret(new CSeqMap(loc));
  598.     if ( scope ) {
  599.         CSeqMap::const_iterator i(
  600.             ret->BeginResolved(scope, size_t(-1), fFindData));
  601.         for ( ; i; ++i ) {
  602.             _ASSERT(i.GetType() == eSeqData);
  603.             switch ( i.GetRefData().Which() ) {
  604.             case CSeq_data::e_Ncbi2na:
  605.             case CSeq_data::e_Ncbi4na:
  606.             case CSeq_data::e_Ncbi8na:
  607.             case CSeq_data::e_Ncbipna:
  608.             case CSeq_data::e_Iupacna:
  609.                 const_cast<CSeqMap&>(*ret).m_Mol = CSeq_inst::eMol_na;
  610.                 break;
  611.             case CSeq_data::e_Ncbi8aa:
  612.             case CSeq_data::e_Ncbieaa:
  613.             case CSeq_data::e_Ncbipaa:
  614.             case CSeq_data::e_Ncbistdaa:
  615.             case CSeq_data::e_Iupacaa:
  616.                 const_cast<CSeqMap&>(*ret).m_Mol = CSeq_inst::eMol_aa;
  617.                 break;
  618.             default:
  619.                 const_cast<CSeqMap&>(*ret).m_Mol = CSeq_inst::eMol_not_set;
  620.             }
  621.         }
  622.     }
  623.     return ret;
  624. }
  627. CConstRef<CSeqMap> CSeqMap::CreateSeqMapForStrand(CConstRef<CSeqMap> seqMap,
  628.                                                   ENa_strand strand)
  629. {
  630.     if ( IsReverse(strand) ) {
  631.         CRef<CSeqMap> ret(new CSeqMap());
  632.         ret->x_AddEnd();
  633.         ret->x_AddSegment(eSeqSubMap,
  634.                           const_cast<CSeqMap*>(seqMap.GetPointer()),
  635.                           0, kInvalidSeqPos, strand);
  636.         ret->x_AddEnd();
  637.         ret->m_Mol = seqMap->m_Mol;
  638.         ret->m_SeqLength = seqMap->m_SeqLength;
  639.         seqMap = ret;
  640.     }
  641.     return seqMap;
  642. }
  645. CSeqMap::CSegment& CSeqMap::x_AddSegment(ESegmentType type, TSeqPos len)
  646. {
  647.     m_Segments.push_back(CSegment(type, len));
  648.     return m_Segments.back();
  649. }
  652. CSeqMap::CSegment& CSeqMap::x_AddSegment(ESegmentType type, TSeqPos len,
  653.                                          const CObject* object)
  654. {
  655.     CSegment& ret = x_AddSegment(type, len);
  656.     ret.m_RefObject.Reset(object);
  657.     return ret;
  658. }
  661. CSeqMap::CSegment& CSeqMap::x_AddSegment(ESegmentType type,
  662.                                          const CObject* object,
  663.                                          TSeqPos refPos,
  664.                                          TSeqPos len,
  665.                                          ENa_strand strand)
  666. {
  667.     CSegment& ret = x_AddSegment(type, len, object);
  668.     ret.m_RefPosition = refPos;
  669.     ret.m_RefMinusStrand = IsReverse(strand);
  670.     return ret;
  671. }
  674. void CSeqMap::x_AddEnd(void)
  675. {
  676.     TSeqPos pos = m_Segments.empty()? 0: kInvalidSeqPos;
  677.     x_AddSegment(eSeqEnd, 0).m_Position = pos;
  678. }
  681. CSeqMap::CSegment& CSeqMap::x_AddGap(TSeqPos len)
  682. {
  683.     return x_AddSegment(eSeqGap, len);
  684. }
  687. CSeqMap::CSegment& CSeqMap::x_AddUnloadedSubMap(TSeqPos len)
  688. {
  689.     return x_AddSegment(eSeqSubMap, len);
  690. }
  693. CSeqMap::CSegment& CSeqMap::x_AddUnloadedSeq_data(TSeqPos len)
  694. {
  695.     return x_AddSegment(eSeqData, len);
  696. }
  699. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_data& data, TSeqPos len)
  700. {
  701.     return x_AddSegment(eSeqData, len, &data);
  702. }
  705. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_point& ref)
  706. {
  707.     return x_AddSegment(eSeqRef, &ref.GetId(),
  708.                         ref.GetPoint(), 1,
  709.                         ref.IsSetStrand()? ref.GetStrand(): eNa_strand_unknown);
  710. }
  713. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_interval& ref)
  714. {
  715.     return x_AddSegment(eSeqRef, &ref.GetId(),
  716.                         ref.GetFrom(), ref.GetLength(),
  717.                         ref.IsSetStrand()? ref.GetStrand(): eNa_strand_unknown);
  718. }
  721. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_id& ref)
  722. {
  723.     return x_AddSegment(eSeqRef, &ref, 0, kInvalidSeqPos);
  724. }
  727. CSeqMap::CSegment& CSeqMap::x_Add(CSeqMap* submap)
  728. {
  729.     return x_AddSegment(eSeqSubMap, kInvalidSeqPos, submap);
  730. }
  733. CSeqMap::CSegment& CSeqMap::x_Add(const CPacked_seqint& seq)
  734. {
  735.     return x_Add(new CSeqMap_Seq_intervals(seq));
  736. }
  739. CSeqMap::CSegment& CSeqMap::x_Add(const CPacked_seqpnt& seq)
  740. {
  741.     return x_Add(new CSeqMap_SeqPoss(seq));
  742. }
  745. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_loc_mix& seq)
  746. {
  747.     return x_Add(new CSeqMap_Seq_locs(seq, seq.Get()));
  748. }
  751. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_loc_equiv& seq)
  752. {
  753.     return x_Add(new CSeqMap_Seq_locs(seq, seq.Get()));
  754. }
  757. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_literal& seq)
  758. {
  759.     if ( seq.IsSetSeq_data() ) {
  760.         return x_Add(seq.GetSeq_data(), seq.GetLength());
  761.     }
  762.     else {
  763.         // No data exist - treat it like a gap
  764.         return x_AddGap(seq.GetLength()); //???
  765.     }
  766. }
  769. CSeqMap::CSegment& CSeqMap::x_Add(const CSeq_loc& loc)
  770. {
  771.     switch ( loc.Which() ) {
  772.     case CSeq_loc::e_not_set:
  773.     case CSeq_loc::e_Null:
  774.     case CSeq_loc::e_Empty:
  775.         return x_AddGap(0); // Add gap ???
  776.     case CSeq_loc::e_Whole:
  777.         return x_Add(loc.GetWhole());
  778.     case CSeq_loc::e_Int:
  779.         return x_Add(loc.GetInt());
  780.     case CSeq_loc::e_Pnt:
  781.         return x_Add(loc.GetPnt());
  782.     case CSeq_loc::e_Packed_int:
  783.         return x_Add(loc.GetPacked_int());
  784.     case CSeq_loc::e_Packed_pnt:
  785.         return x_Add(loc.GetPacked_pnt());
  786.     case CSeq_loc::e_Mix:
  787.         return x_Add(loc.GetMix());
  788.     case CSeq_loc::e_Equiv:
  789.         return x_Add(loc.GetEquiv());
  790.     case CSeq_loc::e_Bond:
  791.         NCBI_THROW(CSeqMapException, eDataError,
  792.                    "e_Bond is not allowed as a reference type");
  793.     case CSeq_loc::e_Feat:
  794.         NCBI_THROW(CSeqMapException, eDataError,
  795.                    "e_Feat is not allowed as a reference type");
  796.     default:
  797.         NCBI_THROW(CSeqMapException, eDataError,
  798.                    "invalid reference type");
  799.     }
  800. }
  803. CSeqMap::CSegment& CSeqMap::x_Add(const CDelta_seq& seq)
  804. {
  805.     switch ( seq.Which() ) {
  806.     case CDelta_seq::e_Loc:
  807.         return x_Add(seq.GetLoc());
  808.     case CDelta_seq::e_Literal:
  809.         return x_Add(seq.GetLiteral());
  810.     default:
  811.         NCBI_THROW(CSeqMapException, eDataError,
  812.                    "Can not add empty Delta-seq");
  813.     }
  814. }
  817. END_SCOPE(objects)
  818. END_NCBI_SCOPE
  820. /*
  821. * ---------------------------------------------------------------------------
  822. * $Log: seq_map.cpp,v $
  823. * Revision 1000.4  2004/06/01 19:24:15  gouriano
  824. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.53
  825. *
  826. * Revision 1.53  2004/05/21 21:42:13  gorelenk
  827. * Added PCH ncbi_pch.hpp
  828. *
  829. * Revision 1.52  2004/03/16 15:47:28  vasilche
  830. * Added CBioseq_set_Handle and set of EditHandles
  831. *
  832. * Revision 1.51  2004/02/19 17:19:10  vasilche
  833. * Removed 'unused argument' warning.
  834. *
  835. * Revision 1.50  2003/11/19 22:18:04  grichenk
  836. * All exceptions are now CException-derived. Catch "exception" rather
  837. * than "runtime_error".
  838. *
  839. * Revision 1.49  2003/11/12 16:53:17  grichenk
  840. * Modified CSeqMap to work with const objects (CBioseq, CSeq_loc etc.)
  841. *
  842. * Revision 1.48  2003/09/30 16:22:04  vasilche
  843. * Updated internal object manager classes to be able to load ID2 data.
  844. * SNP blobs are loaded as ID2 split blobs - readers convert them automatically.
  845. * Scope caches results of requests for data to data loaders.
  846. * Optimized CSeq_id_Handle for gis.
  847. * Optimized bioseq lookup in scope.
  848. * Reduced object allocations in annotation iterators.
  849. * CScope is allowed to be destroyed before other objects using this scope are
  850. * deleted (feature iterators, bioseq handles etc).
  851. * Optimized lookup for matching Seq-ids in CSeq_id_Mapper.
  852. * Added 'adaptive' option to objmgr_demo application.
  853. *
  854. * Revision 1.47  2003/09/05 17:29:40  grichenk
  855. * Structurized Object Manager exceptions
  856. *
  857. * Revision 1.46  2003/08/27 14:27:19  vasilche
  858. * Use Reverse(ENa_strand) function.
  859. *
  860. * Revision 1.45  2003/07/17 19:10:28  grichenk
  861. * Added methods for seq-map and seq-vector validation,
  862. * updated demo.
  863. *
  864. * Revision 1.44  2003/07/14 21:13:26  grichenk
  865. * Added possibility to resolve seq-map iterator withing a single TSE
  866. * and to skip intermediate references during this resolving.
  867. *
  868. * Revision 1.43  2003/06/30 18:39:18  vasilche
  869. * Fixed access to uninitialized member.
  870. *
  871. * Revision 1.42  2003/06/27 19:09:02  grichenk
  872. * Fixed problem with unset sequence length.
  873. *
  874. * Revision 1.41  2003/06/26 19:47:27  grichenk
  875. * Added sequence length cache
  876. *
  877. * Revision 1.40  2003/06/24 14:22:46  vasilche
  878. * Fixed CSeqMap constructor from CSeq_loc.
  879. *
  880. * Revision 1.39  2003/06/12 17:04:55  vasilche
  881. * Fixed creation of CSeqMap for sequences with repr == not_set.
  882. *
  883. * Revision 1.38  2003/06/11 19:32:55  grichenk
  884. * Added molecule type caching to CSeqMap, simplified
  885. * coding and sequence type calculations in CSeqVector.
  886. *
  887. * Revision 1.37  2003/06/02 16:06:38  dicuccio
  888. * Rearranged src/objects/ subtree.  This includes the following shifts:
  889. *     - src/objects/asn2asn --> arc/app/asn2asn
  890. *     - src/objects/testmedline --> src/objects/ncbimime/test
  891. *     - src/objects/objmgr --> src/objmgr
  892. *     - src/objects/util --> src/objmgr/util
  893. *     - src/objects/alnmgr --> src/objtools/alnmgr
  894. *     - src/objects/flat --> src/objtools/flat
  895. *     - src/objects/validator --> src/objtools/validator
  896. *     - src/objects/cddalignview --> src/objtools/cddalignview
  897. * In addition, libseq now includes six of the objects/seq... libs, and libmmdb
  898. * replaces the three libmmdb? libs.
  899. *
  900. * Revision 1.36  2003/05/21 16:03:08  vasilche
  901. * Fixed access to uninitialized optional members.
  902. * Added initialization of mandatory members.
  903. *
  904. * Revision 1.35  2003/05/20 20:36:14  vasilche
  905. * Added FindResolved() with strand argument.
  906. *
  907. * Revision 1.34  2003/05/20 15:44:38  vasilche
  908. * Fixed interaction of CDataSource and CDataLoader in multithreaded app.
  909. * Fixed some warnings on WorkShop.
  910. * Added workaround for memory leak on WorkShop.
  911. *
  912. * Revision 1.33  2003/04/24 16:12:38  vasilche
  913. * Object manager internal structures are splitted more straightforward.
  914. * Removed excessive header dependencies.
  915. *
  916. * Revision 1.32  2003/02/24 18:57:22  vasilche
  917. * Make feature gathering in one linear pass using CSeqMap iterator.
  918. * Do not use feture index by sub locations.
  919. * Sort features at the end of gathering in one vector.
  920. * Extracted some internal structures and classes in separate header.
  921. * Delay creation of mapped features.
  922. *
  923. * Revision 1.31  2003/02/05 17:59:17  dicuccio
  924. * Moved formerly private headers into include/objects/objmgr/impl
  925. *
  926. * Revision 1.30  2003/02/05 15:55:26  vasilche
  927. * Added eSeqEnd segment at the beginning of seq map.
  928. * Added flags to CSeqMap_CI to stop on data, gap, or references.
  929. *
  930. * Revision 1.29  2003/01/28 17:16:06  vasilche
  931. * Added CSeqMap::ResolvedRangeIterator with strand coordinate translation.
  932. *
  933. * Revision 1.28  2003/01/22 20:11:54  vasilche
  934. * Merged functionality of CSeqMapResolved_CI to CSeqMap_CI.
  935. * CSeqMap_CI now supports resolution and iteration over sequence range.
  936. * Added several caches to CScope.
  937. * Optimized CSeqVector().
  938. * Added serveral variants of CBioseqHandle::GetSeqVector().
  939. * Tried to optimize annotations iterator (not much success).
  940. * Rewritten CHandleRange and CHandleRangeMap classes to avoid sorting of list.
  941. *
  942. * Revision 1.27  2002/12/26 20:55:18  dicuccio
  943. * Moved seq_id_mapper.hpp, tse_info.hpp, and bioseq_info.hpp -> include/ tree
  944. *
  945. * Revision 1.26  2002/12/26 20:35:14  ucko
  946. * #include <algorithm> for upper_bound<>
  947. *
  948. * Revision 1.25  2002/12/26 16:39:24  vasilche
  949. * Object manager class CSeqMap rewritten.
  950. *
  951. * Revision 1.24  2002/11/04 21:29:12  grichenk
  952. * Fixed usage of const CRef<> and CRef<> constructor
  953. *
  954. * Revision 1.23  2002/10/18 19:12:40  grichenk
  955. * Removed mutex pools, converted most static mutexes to non-static.
  956. * Protected CSeqMap::x_Resolve() with mutex. Modified code to prevent
  957. * dead-locks.
  958. *
  959. * Revision 1.22  2002/07/08 20:51:02  grichenk
  960. * Moved log to the end of file
  961. * Replaced static mutex (in CScope, CDataSource) with the mutex
  962. * pool. Redesigned CDataSource data locking.
  963. *
  964. * Revision 1.21  2002/05/29 21:21:13  gouriano
  965. * added debug dump
  966. *
  967. * Revision 1.20  2002/05/06 17:43:06  ivanov
  968. * ssize_t changed to long
  969. *
  970. * Revision 1.19  2002/05/06 17:03:49  ivanov
  971. * Sorry. Rollback to R1.17
  972. *
  973. * Revision 1.18  2002/05/06 16:56:23  ivanov
  974. * Fixed typo ssize_t -> size_t
  975. *
  976. * Revision 1.17  2002/05/06 03:28:47  vakatov
  977. * OM/OM1 renaming
  978. *
  979. * Revision 1.16  2002/05/03 21:28:10  ucko
  980. * Introduce T(Signed)SeqPos.
  981. *
  982. * Revision 1.15  2002/05/02 20:42:38  grichenk
  983. * throw -> THROW1_TRACE
  984. *
  985. * Revision 1.14  2002/04/30 18:55:41  gouriano
  986. * added GetRefSeqid function
  987. *
  988. * Revision 1.13  2002/04/11 12:07:30  grichenk
  989. * Redesigned CAnnotTypes_CI to resolve segmented sequences correctly.
  990. *
  991. * Revision 1.12  2002/04/04 21:33:55  grichenk
  992. * Fixed x_FindSegment() for sequences with unresolved segments
  993. *
  994. * Revision 1.11  2002/04/03 18:06:48  grichenk
  995. * Fixed segmented sequence bugs (invalid positioning of literals
  996. * and gaps). Improved CSeqVector performance.
  997. *
  998. * Revision 1.9  2002/03/08 21:36:21  gouriano
  999. * *** empty log message ***
  1000. *
  1001. * Revision 1.8  2002/03/08 21:24:35  gouriano
  1002. * fixed errors with unresolvable references
  1003. *
  1004. * Revision 1.7  2002/02/25 21:05:29  grichenk
  1005. * Removed seq-data references caching. Increased MT-safety. Fixed typos.
  1006. *
  1007. * Revision 1.6  2002/02/21 19:27:06  grichenk
  1008. * Rearranged includes. Added scope history. Added searching for the
  1009. * best seq-id match in data sources and scopes. Updated tests.
  1010. *
  1011. * Revision 1.5  2002/02/01 21:49:38  gouriano
  1012. * minor changes to make it compilable and run on Solaris Workshop
  1013. *
  1014. * Revision 1.4  2002/01/30 22:09:28  gouriano
  1015. * changed CSeqMap interface
  1016. *
  1017. * Revision 1.3  2002/01/23 21:59:32  grichenk
  1018. * Redesigned seq-id handles and mapper
  1019. *
  1020. * Revision 1.2  2002/01/16 16:25:56  gouriano
  1021. * restructured objmgr
  1022. *
  1023. * Revision 1.1  2002/01/11 19:06:24  gouriano
  1024. * restructured objmgr
  1025. *
  1026. *
  1027. * ===========================================================================
  1028. */