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

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

生物技术

开发平台：

C/C++

seq_loc_mapper.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: seq_loc_mapper.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 19:24:12  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: seq_loc_mapper.cpp,v 1000.1 2004/06/01 19:24:12 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. * Author: Aleksey Grichenko
  36. *
  37. * File Description:
  38. *   Seq-loc mapper
  39. *
  40. */
  42. #include <ncbi_pch.hpp>
  43. #include <objmgr/seq_loc_mapper.hpp>
  44. #include <objmgr/scope.hpp>
  45. #include <objmgr/objmgr_exception.hpp>
  46. #include <objmgr/seq_map.hpp>
  47. #include <objmgr/seq_map_ci.hpp>
  48. #include <objmgr/impl/synonyms.hpp>
  49. #include <objmgr/impl/seq_align_mapper.hpp>
  50. #include <objmgr/impl/seq_loc_cvt.hpp>
  51. #include <objects/seqloc/Seq_loc.hpp>
  52. #include <objects/seqfeat/Seq_feat.hpp>
  53. #include <objects/seqfeat/Cdregion.hpp>
  54. #include <objects/seqloc/Seq_loc_equiv.hpp>
  55. #include <objects/seqloc/Seq_bond.hpp>
  56. #include <objects/seqalign/seqalign__.hpp>
  57. #include <algorithm>
  60. BEGIN_NCBI_SCOPE
  61. BEGIN_SCOPE(objects)
  64. CMappingRange::CMappingRange(CSeq_id_Handle     src_id,
  65.                              TSeqPos            src_from,
  66.                              TSeqPos            src_length,
  67.                              ENa_strand         src_strand,
  68.                              CSeq_id_Handle     dst_id,
  69.                              TSeqPos            dst_from,
  70.                              ENa_strand         dst_strand)
  71.     : m_Src_id_Handle(src_id),
  72.       m_Src_from(src_from),
  73.       m_Src_to(src_from + src_length - 1),
  74.       m_Src_strand(src_strand),
  75.       m_Dst_id_Handle(dst_id),
  76.       m_Dst_from(dst_from),
  77.       m_Dst_strand(dst_strand),
  78.       m_Reverse(!SameOrientation(src_strand, dst_strand))
  79. {
  80.     return;
  81. }
  84. bool CMappingRange::CanMap(TSeqPos from,
  85.                            TSeqPos to,
  86.                            bool is_set_strand,
  87.                            ENa_strand strand) const
  88. {
  89.     if (from > m_Src_to || to < m_Src_from) {
  90.         return false;
  91.     }
  92.     return SameOrientation(is_set_strand ? strand : eNa_strand_unknown,
  93.         m_Src_strand)  ||  (m_Src_strand == eNa_strand_unknown);
  94. }
  97. TSeqPos CMappingRange::Map_Pos(TSeqPos pos) const
  98. {
  99.     _ASSERT(pos >= m_Src_from  &&  pos <= m_Src_to);
  100.     if (!m_Reverse) {
  101.         return m_Dst_from + pos - m_Src_from;
  102.     }
  103.     else {
  104.         return m_Dst_from + m_Src_to - pos;
  105.     }
  106. }
  109. CMappingRange::TRange CMappingRange::Map_Range(TSeqPos from, TSeqPos to) const
  110. {
  111.     if (!m_Reverse) {
  112.         return TRange(Map_Pos(max(from, m_Src_from)),
  113.             Map_Pos(min(to, m_Src_to)));
  114.     }
  115.     else {
  116.         return TRange(Map_Pos(min(to, m_Src_to)),
  117.             Map_Pos(max(from, m_Src_from)));
  118.     }
  119. }
  122. bool CMappingRange::Map_Strand(bool is_set_strand,
  123.                                ENa_strand src,
  124.                                ENa_strand* dst) const
  125. {
  126.     _ASSERT(dst);
  127.     if (is_set_strand) {
  128.         *dst = m_Reverse ? Reverse(src) : src;
  129.         return true;
  130.     }
  131.     if (m_Dst_strand != eNa_strand_unknown) {
  132.         // Destination strand may be set for nucleotides
  133.         // even if the source one is not set.
  134.         *dst = m_Dst_strand;
  135.         return true;
  136.     }
  137.     return false;
  138. }
  141. const CMappingRange::TFuzz kEmptyFuzz(0);
  143. CInt_fuzz::ELim CMappingRange::x_ReverseFuzzLim(CInt_fuzz::ELim lim) const
  144. {
  145.     switch ( lim ) {
  146.     case CInt_fuzz::eLim_gt:
  147.         return CInt_fuzz::eLim_lt;
  148.     case CInt_fuzz::eLim_lt:
  149.         return CInt_fuzz::eLim_gt;
  150.     case CInt_fuzz::eLim_tr:
  151.         return CInt_fuzz::eLim_tl;
  152.     case CInt_fuzz::eLim_tl:
  153.         return CInt_fuzz::eLim_tr;
  154.     default:
  155.         return lim;
  156.     }
  157. }
  160. CMappingRange::TRangeFuzz CMappingRange::Map_Fuzz(TRangeFuzz& fuzz) const
  161. {
  162.     // Maps some fuzz types to reverse strand
  163.     if ( !m_Reverse ) {
  164.         return fuzz;
  165.     }
  166.     TRangeFuzz res(fuzz.second, fuzz.first);
  167.     if ( res.first ) {
  168.         switch ( res.first->Which() ) {
  169.         case CInt_fuzz::e_Lim:
  170.             {
  171.                 res.first->SetLim(x_ReverseFuzzLim(res.first->GetLim()));
  172.                 break;
  173.             }
  174.         default:
  175.             // Other types are not converted
  176.             break;
  177.         }
  178.     }
  179.     if ( fuzz.second ) {
  180.         switch ( res.second->Which() ) {
  181.         case CInt_fuzz::e_Lim:
  182.             {
  183.                 res.second->SetLim(x_ReverseFuzzLim(res.second->GetLim()));
  184.                 break;
  185.             }
  186.         default:
  187.             // Other types are not converted
  188.             break;
  189.         }
  190.     }
  191.     return res;
  192. }
  195. CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_feat&  map_feat,
  196.                                  EFeatMapDirection dir,
  197.                                  CScope*           scope)
  198.     : m_Scope(scope),
  199.       m_MergeFlag(eMergeNone),
  200.       m_GapFlag(eGapPreserve),
  201.       m_KeepNonmapping(false),
  202.       m_UseWidth(false),
  203.       m_Dst_width(0)
  204. {
  205.     _ASSERT(map_feat.IsSetProduct());
  206.     int frame = 0;
  207.     if (map_feat.GetData().IsCdregion()) {
  208.         frame = map_feat.GetData().GetCdregion().GetFrame();
  209.     }
  210.     if (dir == eLocationToProduct) {
  211.         x_Initialize(map_feat.GetLocation(), map_feat.GetProduct(), frame);
  212.     }
  213.     else {
  214.         x_Initialize(map_feat.GetProduct(), map_feat.GetLocation(), frame);
  215.     }
  216. }
  219. CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_loc& source,
  220.                                  const CSeq_loc& target,
  221.                                  CScope* scope)
  222.     : m_Scope(scope),
  223.       m_MergeFlag(eMergeNone),
  224.       m_GapFlag(eGapPreserve),
  225.       m_KeepNonmapping(false),
  226.       m_UseWidth(false),
  227.       m_Dst_width(0)
  228. {
  229.     x_Initialize(source, target);
  230. }
  233. CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_align& map_align,
  234.                                  const CSeq_id&    to_id,
  235.                                  CScope*           scope)
  236.     : m_Scope(scope),
  237.       m_MergeFlag(eMergeNone),
  238.       m_GapFlag(eGapPreserve),
  239.       m_KeepNonmapping(false),
  240.       m_UseWidth(false),
  241.       m_Dst_width(0)
  242. {
  243.     x_Initialize(map_align, to_id);
  244. }
  247. CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_align& map_align,
  248.                                  size_t            to_row,
  249.                                  CScope*           scope)
  250.     : m_Scope(scope),
  251.       m_MergeFlag(eMergeNone),
  252.       m_GapFlag(eGapPreserve),
  253.       m_KeepNonmapping(false),
  254.       m_UseWidth(false),
  255.       m_Dst_width(0)
  256. {
  257.     x_Initialize(map_align, to_row);
  258. }
  261. CSeq_loc_Mapper::CSeq_loc_Mapper(CBioseq_Handle target_seq)
  262.     : m_Scope(&target_seq.GetScope()),
  263.       m_MergeFlag(eMergeNone),
  264.       m_GapFlag(eGapPreserve),
  265.       m_KeepNonmapping(false),
  266.       m_UseWidth(false),
  267.       m_Dst_width(0)
  268. {
  269.     CConstRef<CSeq_id> dst_id = target_seq.GetSeqId();
  270.     x_Initialize(target_seq.GetSeqMap(),
  271.         dst_id.GetPointerOrNull());
  272.     // Ignore seq-map destination ranges, map whole sequence to itself,
  273.     // use unknown strand only.
  274.     m_DstRanges.resize(1);
  275.     m_DstRanges[0].clear();
  276.     if (m_Scope) {
  277.         CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*dst_id);
  278.         ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
  279.             m_DstRanges[0][target_seq.GetSeq_id_Handle()]
  280.                 .push_back(TRange::GetWhole());
  281.         }
  282.     }
  283.     else {
  284.         m_DstRanges[0][target_seq.GetSeq_id_Handle()]
  285.             .push_back(TRange::GetWhole());
  286.     }
  287. }
  290. CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeqMap& seq_map,
  291.                                  const CSeq_id* dst_id,
  292.                                  CScope*        scope)
  293.     : m_Scope(scope),
  294.       m_MergeFlag(eMergeNone),
  295.       m_GapFlag(eGapPreserve),
  296.       m_KeepNonmapping(false),
  297.       m_UseWidth(false),
  298.       m_Dst_width(0)
  299. {
  300.     x_Initialize(seq_map, dst_id);
  301. }
  304. CSeq_loc_Mapper::CSeq_loc_Mapper(size_t          depth,
  305.                                  CBioseq_Handle& source_seq)
  306.     : m_Scope(&source_seq.GetScope()),
  307.       m_MergeFlag(eMergeNone),
  308.       m_GapFlag(eGapPreserve),
  309.       m_KeepNonmapping(false),
  310.       m_UseWidth(false),
  311.       m_Dst_width(0)
  312. {
  313.     if (depth > 0) {
  314.         depth--;
  315.         x_Initialize(source_seq.GetSeqMap(),
  316.             depth, source_seq.GetSeqId().GetPointer());
  317.     }
  318.     else /* if (depth == 0) */ {
  319.         // Synonyms conversion
  320.         CConstRef<CSeq_id> dst_id = source_seq.GetSeqId();
  321.         m_DstRanges.resize(1);
  322.         if (m_Scope) {
  323.             CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*dst_id);
  324.             ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
  325.                 m_DstRanges[0][source_seq.GetSeq_id_Handle()]
  326.                     .push_back(TRange::GetWhole());
  327.             }
  328.         }
  329.         else {
  330.             m_DstRanges[0][source_seq.GetSeq_id_Handle()]
  331.                 .push_back(TRange::GetWhole());
  332.         }
  333.     }
  334. }
  337. CSeq_loc_Mapper::CSeq_loc_Mapper(size_t         depth,
  338.                                  const CSeqMap& source_seqmap,
  339.                                  const CSeq_id* src_id,
  340.                                  CScope*        scope)
  341.     : m_Scope(scope),
  342.       m_MergeFlag(eMergeNone),
  343.       m_GapFlag(eGapPreserve),
  344.       m_KeepNonmapping(false),
  345.       m_UseWidth(false),
  346.       m_Dst_width(0)
  347. {
  348.     if (depth > 0) {
  349.         depth--;
  350.         x_Initialize(source_seqmap, depth, src_id);
  351.     }
  352.     else /* if (depth == 0) */ {
  353.         // Synonyms conversion
  354.         m_DstRanges.resize(1);
  355.         if (bool(m_Scope)  &&  src_id) {
  356.             CSeq_id_Handle src_idh = CSeq_id_Handle::GetHandle(*src_id);
  357.             CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*src_id);
  358.             ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
  359.                 m_DstRanges[0][src_idh]
  360.                     .push_back(TRange::GetWhole());
  361.             }
  362.         }
  363.     }
  364. }
  367. CSeq_loc_Mapper::~CSeq_loc_Mapper(void)
  368. {
  369.     return;
  370. }
  373. void CSeq_loc_Mapper::PreserveDestinationLocs(void)
  374. {
  375.     for (size_t str_idx = 0; str_idx < m_DstRanges.size(); str_idx++) {
  376.         NON_CONST_ITERATE(TDstIdMap, id_it, m_DstRanges[str_idx]) {
  377.             id_it->second.sort();
  378.             TSeqPos dst_start = kInvalidSeqPos;
  379.             TSeqPos dst_stop = kInvalidSeqPos;
  380.             ITERATE(TDstRanges, rg_it, id_it->second) {
  381.                 // Collect and merge ranges
  382.                 if (dst_start == kInvalidSeqPos) {
  383.                     dst_start = rg_it->GetFrom();
  384.                     dst_stop = rg_it->GetTo();
  385.                     continue;
  386.                 }
  387.                 if (rg_it->GetFrom() <= dst_stop + 1) {
  388.                     // overlapping or abutting ranges, continue collecting
  389.                     dst_stop = max(dst_stop, rg_it->GetTo());
  390.                     continue;
  391.                 }
  392.                 // Separate ranges, add conversion and restart collecting
  393.                 CRef<CMappingRange> cvt(new CMappingRange(
  394.                     id_it->first, dst_start, dst_stop - dst_start + 1,
  395.                     ENa_strand(str_idx),
  396.                     id_it->first, dst_start, ENa_strand(str_idx)));
  397.                 m_IdMap[cvt->m_Src_id_Handle].insert(TRangeMap::value_type(
  398.                     TRange(cvt->m_Src_from, cvt->m_Src_to), cvt));
  399.             }
  400.             if (dst_start < dst_stop) {
  401.                 CRef<CMappingRange> cvt(new CMappingRange(
  402.                     id_it->first, dst_start, dst_stop - dst_start + 1,
  403.                     ENa_strand(str_idx),
  404.                     id_it->first, dst_start, ENa_strand(str_idx)));
  405.                 m_IdMap[cvt->m_Src_id_Handle].insert(TRangeMap::value_type(
  406.                     TRange(cvt->m_Src_from, cvt->m_Src_to), cvt));
  407.             }
  408.         }
  409.     }
  410.     m_DstRanges.clear();
  411. }
  414. CSeq_loc_Mapper::TMappedRanges&
  415. CSeq_loc_Mapper::x_GetMappedRanges(const CSeq_id_Handle& id,
  416.                                    int strand_idx) const
  417. {
  418.     TRangesByStrand& str_vec = m_MappedLocs[id];
  419.     if ((int)str_vec.size() <= strand_idx) {
  420.         str_vec.resize(strand_idx + 1);
  421.     }
  422.     return str_vec[strand_idx];
  423. }
  426. void CSeq_loc_Mapper::x_PushRangesToDstMix(void)
  427. {
  428.     if (m_MappedLocs.size() == 0) {
  429.         return;
  430.     }
  431.     // Push everything already mapped to the destination mix
  432.     CRef<CSeq_loc> loc = x_GetMappedSeq_loc();
  433.     if ( !m_Dst_loc ) {
  434.         m_Dst_loc = loc;
  435.         return;
  436.     }
  437.     if ( !loc->IsNull() ) {
  438.         x_PushLocToDstMix(loc);
  439.     }
  440. }
  443. int CSeq_loc_Mapper::x_CheckSeqWidth(const CSeq_id& id, int width)
  444. {
  445.     if (!m_Scope) {
  446.         return width;
  447.     }
  448.     CBioseq_Handle handle;
  449.     try {
  450.         handle = m_Scope->GetBioseqHandle(id);
  451.     } catch (...) {
  452.         return width;
  453.     }
  454.     if ( !handle ) {
  455.         return width;
  456.     }
  457.     switch ( handle.GetBioseqMolType() ) {
  458.     case CSeq_inst::eMol_dna:
  459.     case CSeq_inst::eMol_rna:
  460.     case CSeq_inst::eMol_na:
  461.         {
  462.             if ( width != 3 ) {
  463.                 if ( width ) {
  464.                     // width already set to a different value
  465.                     NCBI_THROW(CLocMapperException, eBadLocation,
  466.                                "Location contains different sequence types");
  467.                 }
  468.                 width = 3;
  469.             }
  470.             break;
  471.         }
  472.     case CSeq_inst::eMol_aa:
  473.         {
  474.             if ( width != 1 ) {
  475.                 if ( width ) {
  476.                     // width already set to a different value
  477.                     NCBI_THROW(CLocMapperException, eBadLocation,
  478.                                "Location contains different sequence types");
  479.                 }
  480.                 width = 1;
  481.             }
  482.             break;
  483.         }
  484.     default:
  485.         return width;
  486.     }
  487.     // Destination width should always be checked first
  488.     if ( !m_Dst_width ) {
  489.         m_Dst_width = width;
  490.     }
  491.     TWidthFlags wid_cvt = 0;
  492.     if (width == 1) {
  493.         wid_cvt |= fWidthProtToNuc;
  494.     }
  495.     if (m_Dst_width == 1) {
  496.         wid_cvt |= fWidthNucToProt;
  497.     }
  498.     CConstRef<CSynonymsSet> syns = m_Scope->GetSynonyms(id);
  499.     ITERATE(CSynonymsSet, syn_it, *syns) {
  500.         m_Widths[syns->GetSeq_id_Handle(syn_it)] = wid_cvt;
  501.     }
  502.     return width;
  503. }
  506. int CSeq_loc_Mapper::x_CheckSeqWidth(const CSeq_loc& loc,
  507.                                      TSeqPos* total_length)
  508. {
  509.     int width = 0;
  510.     *total_length = 0;
  511.     for (CSeq_loc_CI it(loc); it; ++it) {
  512.         *total_length += it.GetRange().GetLength();
  513.         if ( !m_Scope ) {
  514.             continue; // don't check sequence type;
  515.         }
  516.         width = x_CheckSeqWidth(it.GetSeq_id(), width);
  517.     }
  518.     return width;
  519. }
  522. TSeqPos CSeq_loc_Mapper::x_GetRangeLength(const CSeq_loc_CI& it)
  523. {
  524.     if (it.IsWhole()  &&  IsReverse(it.GetStrand())) {
  525.         // For reverse strand we need real interval length, not just "whole"
  526.         CBioseq_Handle h;
  527.         if (m_Scope) {
  528.             h = m_Scope->GetBioseqHandle(it.GetSeq_id());
  529.         }
  530.         if ( !h ) {
  531.             NCBI_THROW(CLocMapperException, eUnknownLength,
  532.                        "Can not map from minus strand -- "
  533.                        "unknown sequence length");
  534.         }
  535.         return h.GetBioseqLength();
  536.     }
  537.     else {
  538.         return it.GetRange().GetLength();
  539.     }
  540. }
  543. void CSeq_loc_Mapper::x_AddConversion(const CSeq_id& src_id,
  544.                                       TSeqPos        src_start,
  545.                                       ENa_strand     src_strand,
  546.                                       const CSeq_id& dst_id,
  547.                                       TSeqPos        dst_start,
  548.                                       ENa_strand     dst_strand,
  549.                                       TSeqPos        length)
  550. {
  551.     if (m_DstRanges.size() <= size_t(dst_strand)) {
  552.         m_DstRanges.resize(size_t(dst_strand) + 1);
  553.     }
  554.     if (m_Scope) {
  555.         CConstRef<CSynonymsSet> syns = m_Scope->GetSynonyms(src_id);
  556.         ITERATE(CSynonymsSet, syn_it, *syns) {
  557.             CRef<CMappingRange> cvt(new CMappingRange(
  558.                 CSynonymsSet::GetSeq_id_Handle(syn_it),
  559.                 src_start, length, src_strand,
  560.                 CSeq_id_Handle::GetHandle(dst_id), dst_start, dst_strand));
  561.             m_IdMap[cvt->m_Src_id_Handle].insert(TRangeMap::value_type(
  562.                 TRange(cvt->m_Src_from, cvt->m_Src_to), cvt));
  563.         }
  564.         CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(dst_id);
  565.         ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
  566.             m_DstRanges[size_t(dst_strand)]
  567.                 [CSynonymsSet::GetSeq_id_Handle(dst_syn_it)]
  568.                 .push_back(TRange(dst_start, dst_start + length - 1));
  569.         }
  570.     }
  571.     else {
  572.         CRef<CMappingRange> cvt(new CMappingRange(
  573.             CSeq_id_Handle::GetHandle(src_id),
  574.             src_start, length, src_strand,
  575.             CSeq_id_Handle::GetHandle(dst_id), dst_start, dst_strand));
  576.         m_IdMap[cvt->m_Src_id_Handle].insert(TRangeMap::value_type(
  577.             TRange(cvt->m_Src_from, cvt->m_Src_to), cvt));
  578.         m_DstRanges[size_t(dst_strand)][CSeq_id_Handle::GetHandle(dst_id)]
  579.             .push_back(TRange(dst_start, dst_start + length - 1));
  580.     }
  581. }
  584. void CSeq_loc_Mapper::x_NextMappingRange(const CSeq_id& src_id,
  585.                                          TSeqPos&       src_start,
  586.                                          TSeqPos&       src_len,
  587.                                          ENa_strand     src_strand,
  588.                                          const CSeq_id& dst_id,
  589.                                          TSeqPos&       dst_start,
  590.                                          TSeqPos&       dst_len,
  591.                                          ENa_strand     dst_strand)
  592. {
  593.     TSeqPos cvt_src_start = src_start;
  594.     TSeqPos cvt_dst_start = dst_start;
  595.     TSeqPos cvt_length;
  597.     if (src_len == dst_len) {
  598.         cvt_length = src_len;
  599.         src_len = 0;
  600.         dst_len = 0;
  601.     }
  602.     else if (src_len > dst_len) {
  603.         // reverse interval for minus strand
  604.         if (IsReverse(src_strand)) {
  605.             cvt_src_start += src_len - dst_len;
  606.         }
  607.         else {
  608.             src_start += dst_len;
  609.         }
  610.         cvt_length = dst_len;
  611.         src_len -= cvt_length;
  612.         dst_len = 0;
  613.     }
  614.     else { // if (src_len < dst_len)
  615.         // reverse interval for minus strand
  616.         if ( IsReverse(dst_strand) ) {
  617.             cvt_dst_start += dst_len - src_len;
  618.         }
  619.         else {
  620.             dst_start += src_len;
  621.         }
  622.         cvt_length = src_len;
  623.         dst_len -= cvt_length;
  624.         src_len = 0;
  625.     }
  626.     x_AddConversion(
  627.         src_id, cvt_src_start, src_strand,
  628.         dst_id, cvt_dst_start, dst_strand,
  629.         cvt_length);
  630. }
  633. void CSeq_loc_Mapper::x_Initialize(const CSeq_loc& source,
  634.                                    const CSeq_loc& target,
  635.                                    int             frame)
  636. {
  637.     // Check sequence type or total location length to
  638.     // adjust intervals according to character width.
  639.     TSeqPos dst_total_len;
  640.     x_CheckSeqWidth(target, &dst_total_len);
  641.     TSeqPos src_total_len;
  642.     int src_width = x_CheckSeqWidth(source, &src_total_len);
  644.     if (!src_width  ||  !m_Dst_width) {
  645.         // Try to detect types by lengths
  646.         if (src_total_len == kInvalidSeqPos || dst_total_len== kInvalidSeqPos) {
  647.             NCBI_THROW(CLocMapperException, eBadLocation,
  648.                        "Undefined location length -- "
  649.                        "unable to detect sequence type");
  650.         }
  651.         if (src_total_len == dst_total_len) {
  652.             src_width = 1;
  653.             m_Dst_width = 1;
  654.             _ASSERT(!frame);
  655.         }
  656.         // truncate incomplete left and right codons
  657.         else if (src_total_len/3 == dst_total_len) {
  658.             src_width = 3;
  659.             m_Dst_width = 1;
  660.         }
  661.         else if (dst_total_len/3 == src_total_len) {
  662.             src_width = 1;
  663.             m_Dst_width = 3;
  664.         }
  665.         else {
  666.             NCBI_THROW(CAnnotException, eBadLocation,
  667.                        "Wrong location length -- "
  668.                        "unable to detect sequence type");
  669.         }
  670.     }
  671.     else {
  672.         if (src_width == m_Dst_width) {
  673.             src_width = 1;
  674.             m_Dst_width = 1;
  675.             _ASSERT(!frame);
  676.         }
  677.     }
  678.     m_UseWidth |= (src_width != m_Dst_width);
  680.     // Create conversions
  681.     CSeq_loc_CI src_it(source);
  682.     CSeq_loc_CI dst_it(target);
  683.     TSeqPos src_start = src_it.GetRange().GetFrom()*m_Dst_width;
  684.     TSeqPos src_len = x_GetRangeLength(src_it)*m_Dst_width;
  685.     TSeqPos dst_start = dst_it.GetRange().GetFrom()*src_width;
  686.     TSeqPos dst_len = x_GetRangeLength(dst_it)*src_width;
  687.     if ( frame ) {
  688.         // ignore pre-frame range
  689.         if (src_width == 3) {
  690.             src_start += frame - 1;
  691.         }
  692.         if (m_Dst_width == 3) {
  693.             dst_start += frame - 1;
  694.         }
  695.     }
  696.     while (src_it  &&  dst_it) {
  697.         x_NextMappingRange(
  698.             src_it.GetSeq_id(), src_start, src_len, src_it.GetStrand(),
  699.             dst_it.GetSeq_id(), dst_start, dst_len, dst_it.GetStrand());
  700.         if (src_len == 0  &&  ++src_it) {
  701.             src_start = src_it.GetRange().GetFrom()*m_Dst_width;
  702.             src_len = x_GetRangeLength(src_it)*m_Dst_width;
  703.         }
  704.         if (dst_len == 0  &&  ++dst_it) {
  705.             dst_start = dst_it.GetRange().GetFrom()*src_width;
  706.             dst_len = x_GetRangeLength(dst_it)*src_width;
  707.         }
  708.     }
  709. }
  712. void CSeq_loc_Mapper::x_Initialize(const CSeq_align& map_align,
  713.                                    const CSeq_id&    to_id)
  714. {
  715.     switch ( map_align.GetSegs().Which() ) {
  716.     case CSeq_align::C_Segs::e_Dendiag:
  717.         {
  718.             const TDendiag& diags = map_align.GetSegs().GetDendiag();
  719.             ITERATE(TDendiag, diag_it, diags) {
  720.                 size_t to_row = size_t(-1);
  721.                 for (size_t i = 0; i < (*diag_it)->GetIds().size(); ++i) {
  722.                     // find the first row with required ID
  723.                     // ??? check if there are multiple rows with the same ID?
  724.                     if ( (*diag_it)->GetIds()[i]->Equals(to_id) ) {
  725.                         to_row = i;
  726.                         break;
  727.                     }
  728.                 }
  729.                 if (to_row == size_t(-1)) {
  730.                     NCBI_THROW(CLocMapperException, eBadAlignment,
  731.                                "Target ID not found in the alignment");
  732.                 }
  733.                 x_InitAlign(**diag_it, to_row);
  734.             }
  735.             break;
  736.         }
  737.     case CSeq_align::C_Segs::e_Denseg:
  738.         {
  739.             const CDense_seg& dseg = map_align.GetSegs().GetDenseg();
  740.             size_t to_row = size_t(-1);
  741.             for (size_t i = 0; i < dseg.GetIds().size(); ++i) {
  742.                 if (dseg.GetIds()[i]->Equals(to_id)) {
  743.                     to_row = i;
  744.                     break;
  745.                 }
  746.             }
  747.             if (to_row == size_t(-1)) {
  748.                 NCBI_THROW(CLocMapperException, eBadAlignment,
  749.                            "Target ID not found in the alignment");
  750.             }
  751.             x_InitAlign(dseg, to_row);
  752.             break;
  753.         }
  754.     case CSeq_align::C_Segs::e_Std:
  755.         {
  756.             const TStd& std_segs = map_align.GetSegs().GetStd();
  757.             ITERATE(TStd, std_seg, std_segs) {
  758.                 size_t to_row = size_t(-1);
  759.                 for (size_t i = 0; i < (*std_seg)->GetIds().size(); ++i) {
  760.                     if ((*std_seg)->GetIds()[i]->Equals(to_id)) {
  761.                         to_row = i;
  762.                         break;
  763.                     }
  764.                 }
  765.                 if (to_row == size_t(-1)) {
  766.                     NCBI_THROW(CLocMapperException, eBadAlignment,
  767.                                "Target ID not found in the alignment");
  768.                 }
  769.                 x_InitAlign(**std_seg, to_row);
  770.             }
  771.             break;
  772.         }
  773.     case CSeq_align::C_Segs::e_Packed:
  774.         {
  775.             const CPacked_seg& pseg = map_align.GetSegs().GetPacked();
  776.             size_t to_row = size_t(-1);
  777.             for (size_t i = 0; i < pseg.GetIds().size(); ++i) {
  778.                 if (pseg.GetIds()[i]->Equals(to_id)) {
  779.                     to_row = i;
  780.                     break;
  781.                 }
  782.             }
  783.             if (to_row == size_t(-1)) {
  784.                 NCBI_THROW(CLocMapperException, eBadAlignment,
  785.                            "Target ID not found in the alignment");
  786.             }
  787.             x_InitAlign(pseg, to_row);
  788.             break;
  789.         }
  790.     case CSeq_align::C_Segs::e_Disc:
  791.         {
  792.             const CSeq_align_set& aln_set = map_align.GetSegs().GetDisc();
  793.             ITERATE(CSeq_align_set::Tdata, aln, aln_set.Get()) {
  794.                 x_Initialize(**aln, to_id);
  795.             }
  796.             break;
  797.         }
  798.     default:
  799.         NCBI_THROW(CLocMapperException, eBadAlignment,
  800.                    "Unsupported alignment type");
  801.     }
  802. }
  805. void CSeq_loc_Mapper::x_Initialize(const CSeq_align& map_align,
  806.                                    size_t            to_row)
  807. {
  808.     switch ( map_align.GetSegs().Which() ) {
  809.     case CSeq_align::C_Segs::e_Dendiag:
  810.         {
  811.             const TDendiag& diags = map_align.GetSegs().GetDendiag();
  812.             ITERATE(TDendiag, diag_it, diags) {
  813.                 x_InitAlign(**diag_it, to_row);
  814.             }
  815.             break;
  816.         }
  817.     case CSeq_align::C_Segs::e_Denseg:
  818.         {
  819.             const CDense_seg& dseg = map_align.GetSegs().GetDenseg();
  820.             x_InitAlign(dseg, to_row);
  821.             break;
  822.         }
  823.     case CSeq_align::C_Segs::e_Std:
  824.         {
  825.             const TStd& std_segs = map_align.GetSegs().GetStd();
  826.             ITERATE(TStd, std_seg, std_segs) {
  827.                 x_InitAlign(**std_seg, to_row);
  828.             }
  829.             break;
  830.         }
  831.     case CSeq_align::C_Segs::e_Packed:
  832.         {
  833.             const CPacked_seg& pseg = map_align.GetSegs().GetPacked();
  834.             x_InitAlign(pseg, to_row);
  835.             break;
  836.         }
  837.     case CSeq_align::C_Segs::e_Disc:
  838.         {
  839.             // The same row in each sub-alignment?
  840.             const CSeq_align_set& aln_set = map_align.GetSegs().GetDisc();
  841.             ITERATE(CSeq_align_set::Tdata, aln, aln_set.Get()) {
  842.                 x_Initialize(**aln, to_row);
  843.             }
  844.             break;
  845.         }
  846.     default:
  847.         NCBI_THROW(CLocMapperException, eBadAlignment,
  848.                    "Unsupported alignment type");
  849.     }
  850. }
  853. void CSeq_loc_Mapper::x_InitAlign(const CDense_diag& diag, size_t to_row)
  854. {
  855.     if (diag.GetStarts()[to_row] < 0) {
  856.         // Destination ID is not present in this dense-diag
  857.         return;
  858.     }
  859.     size_t dim = diag.GetDim();
  860.     _ASSERT(to_row < dim);
  861.     if (dim != diag.GetIds().size()) {
  862.         ERR_POST(Warning << "Invalid 'ids' size in dendiag");
  863.         dim = min(dim, diag.GetIds().size());
  864.     }
  865.     if (dim != diag.GetStarts().size()) {
  866.         ERR_POST(Warning << "Invalid 'starts' size in dendiag");
  867.         dim = min(dim, diag.GetStarts().size());
  868.     }
  869.     bool have_strands = diag.IsSetStrands();
  870.     if (have_strands && dim != diag.GetStrands().size()) {
  871.         ERR_POST(Warning << "Invalid 'strands' size in dendiag");
  872.         dim = min(dim, diag.GetStrands().size());
  873.     }
  875.     ENa_strand dst_strand = have_strands ?
  876.         diag.GetStrands()[to_row] : eNa_strand_unknown;
  877.     const CSeq_id& dst_id = *diag.GetIds()[to_row];
  878.     if (!m_Dst_width) {
  879.         m_Dst_width = x_CheckSeqWidth(dst_id, m_Dst_width);
  880.     }
  882.     for (size_t row = 0; row < dim; ++row) {
  883.         if (row == to_row) {
  884.             continue;
  885.         }
  886.         if (diag.GetStarts()[row] < 0) {
  887.             continue;
  888.         }
  889.         // In alignments with multiple sequence types segment length
  890.         // should be multiplied by character width.
  891.         int src_width = 0;
  892.         const CSeq_id& src_id = *diag.GetIds()[row];
  893.         src_width = x_CheckSeqWidth(src_id, src_width);
  894.         m_UseWidth |= (src_width != m_Dst_width);
  895.         int dst_width_rel = (m_UseWidth) ? 1 : m_Dst_width;
  896.         int src_width_rel = (m_UseWidth) ? 1 : src_width;
  897.         TSeqPos src_start = diag.GetStarts()[row]*dst_width_rel;
  898.         TSeqPos src_len = diag.GetLen()*src_width*dst_width_rel;
  899.         TSeqPos dst_start = diag.GetStarts()[to_row]*src_width_rel;
  900.         TSeqPos dst_len = diag.GetLen()*m_Dst_width*src_width_rel;
  901.         ENa_strand src_strand = have_strands ?
  902.             diag.GetStrands()[row] : eNa_strand_unknown;
  903.         x_NextMappingRange(
  904.             src_id, src_start, src_len, src_strand,
  905.             dst_id, dst_start, dst_len, dst_strand);
  906.         _ASSERT(!src_len  &&  !dst_len);
  907.     }
  908. }
  911. void CSeq_loc_Mapper::x_InitAlign(const CDense_seg& denseg, size_t to_row)
  912. {
  913.     size_t dim = denseg.GetDim();
  914.     _ASSERT(to_row < dim);
  916.     size_t numseg = denseg.GetNumseg();
  917.     // claimed dimension may not be accurate :-/
  918.     if (numseg != denseg.GetLens().size()) {
  919.         ERR_POST(Warning << "Invalid 'lens' size in denseg");
  920.         numseg = min(numseg, denseg.GetLens().size());
  921.     }
  922.     if (dim != denseg.GetIds().size()) {
  923.         ERR_POST(Warning << "Invalid 'ids' size in denseg");
  924.         dim = min(dim, denseg.GetIds().size());
  925.     }
  926.     if (dim*numseg != denseg.GetStarts().size()) {
  927.         ERR_POST(Warning << "Invalid 'starts' size in denseg");
  928.         dim = min(dim*numseg, denseg.GetStarts().size()) / numseg;
  929.     }
  930.     bool have_strands = denseg.IsSetStrands();
  931.     if (have_strands && dim*numseg != denseg.GetStrands().size()) {
  932.         ERR_POST(Warning << "Invalid 'strands' size in denseg");
  933.         dim = min(dim*numseg, denseg.GetStrands().size()) / numseg;
  934.     }
  936.     const CSeq_id& dst_id = *denseg.GetIds()[to_row];
  938.     if (!m_Dst_width) {
  939.         if ( denseg.IsSetWidths() ) {
  940.             m_Dst_width = denseg.GetWidths()[to_row];
  941.         }
  942.         else {
  943.             m_Dst_width = x_CheckSeqWidth(dst_id, m_Dst_width);
  944.         }
  945.     }
  946.     for (size_t row = 0; row < dim; ++row) {
  947.         if (row == to_row) {
  948.             continue;
  949.         }
  950.         const CSeq_id& src_id = *denseg.GetIds()[row];
  952.         int src_width = 0;
  953.         if ( denseg.IsSetWidths() ) {
  954.             src_width = denseg.GetWidths()[row];
  955.         }
  956.         else {
  957.             src_width = x_CheckSeqWidth(src_id, src_width);
  958.         }
  959.         m_UseWidth |= (src_width != m_Dst_width);
  960.         int dst_width_rel = (m_UseWidth) ? 1 : m_Dst_width;
  961.         int src_width_rel = (m_UseWidth) ? 1 : src_width;
  963.         for (size_t seg = 0; seg < numseg; ++seg) {
  964.             int i_src_start = denseg.GetStarts()[seg*dim + row];
  965.             int i_dst_start = denseg.GetStarts()[seg*dim + to_row];
  966.             if (i_src_start < 0  ||  i_dst_start < 0) {
  967.                 continue;
  968.             }
  970.             ENa_strand dst_strand = have_strands ?
  971.                 denseg.GetStrands()[seg*dim + to_row] : eNa_strand_unknown;
  972.             ENa_strand src_strand = have_strands ?
  973.                 denseg.GetStrands()[seg*dim + row] : eNa_strand_unknown;
  975.             // In alignments with multiple sequence types segment length
  976.             // should be multiplied by character width.
  977.             TSeqPos src_len = denseg.GetLens()[seg]*src_width*dst_width_rel;
  978.             TSeqPos dst_len = denseg.GetLens()[seg]*m_Dst_width*src_width_rel;
  979.             TSeqPos src_start = (TSeqPos)(i_src_start)*dst_width_rel;
  980.             TSeqPos dst_start = (TSeqPos)(i_dst_start)*src_width_rel;
  981.             x_NextMappingRange(
  982.                 src_id, src_start, src_len, src_strand,
  983.                 dst_id, dst_start, dst_len, dst_strand);
  984.             _ASSERT(!src_len  &&  !dst_len);
  985.         }
  986.     }
  987. }
  990. void CSeq_loc_Mapper::x_InitAlign(const CStd_seg& sseg, size_t to_row)
  991. {
  992.     size_t dim = sseg.GetDim();
  993.     if (dim != sseg.GetLoc().size()) {
  994.         ERR_POST(Warning << "Invalid 'loc' size in std-seg");
  995.         dim = min(dim, sseg.GetLoc().size());
  996.     }
  997.     if (sseg.IsSetIds()
  998.         && dim != sseg.GetIds().size()) {
  999.         ERR_POST(Warning << "Invalid 'ids' size in std-seg");
  1000.         dim = min(dim, sseg.GetIds().size());
  1001.     }
  1003.     const CSeq_loc& dst_loc = *sseg.GetLoc()[to_row];
  1005.     for (size_t row = 0; row < dim; ++row ) {
  1006.         if (row == to_row) {
  1007.             continue;
  1008.         }
  1009.         const CSeq_loc& src_loc = *sseg.GetLoc()[row];
  1010.         if ( src_loc.IsEmpty() ) {
  1011.             // skipped row in this segment
  1012.             continue;
  1013.         }
  1014.         x_Initialize(src_loc, dst_loc);
  1015.     }
  1016. }
  1019. void CSeq_loc_Mapper::x_InitAlign(const CPacked_seg& pseg, size_t to_row)
  1020. {
  1021.     size_t dim    = pseg.GetDim();
  1022.     size_t numseg = pseg.GetNumseg();
  1023.     // claimed dimension may not be accurate :-/
  1024.     if (numseg != pseg.GetLens().size()) {
  1025.         ERR_POST(Warning << "Invalid 'lens' size in packed-seg");
  1026.         numseg = min(numseg, pseg.GetLens().size());
  1027.     }
  1028.     if (dim != pseg.GetIds().size()) {
  1029.         ERR_POST(Warning << "Invalid 'ids' size in packed-seg");
  1030.         dim = min(dim, pseg.GetIds().size());
  1031.     }
  1032.     if (dim*numseg != pseg.GetStarts().size()) {
  1033.         ERR_POST(Warning << "Invalid 'starts' size in packed-seg");
  1034.         dim = min(dim*numseg, pseg.GetStarts().size()) / numseg;
  1035.     }
  1036.     bool have_strands = pseg.IsSetStrands();
  1037.     if (have_strands && dim*numseg != pseg.GetStrands().size()) {
  1038.         ERR_POST(Warning << "Invalid 'strands' size in packed-seg");
  1039.         dim = min(dim*numseg, pseg.GetStrands().size()) / numseg;
  1040.     }
  1042.     const CSeq_id& dst_id = *pseg.GetIds()[to_row];
  1043.     if (!m_Dst_width) {
  1044.         m_Dst_width = x_CheckSeqWidth(dst_id, m_Dst_width);
  1045.     }
  1046.     for (size_t row = 0; row < dim; ++row) {
  1047.         if (row == to_row) {
  1048.             continue;
  1049.         }
  1050.         const CSeq_id& src_id = *pseg.GetIds()[row];
  1052.         int src_width = 0;
  1053.         src_width = x_CheckSeqWidth(src_id, src_width);
  1054.         m_UseWidth |= (src_width != m_Dst_width);
  1055.         int dst_width_rel = (m_UseWidth) ? 1 : m_Dst_width;
  1056.         int src_width_rel = (m_UseWidth) ? 1 : src_width;
  1058.         for (size_t seg = 0; seg < numseg; ++seg) {
  1059.             if (!pseg.GetPresent()[row]  ||  !pseg.GetPresent()[to_row]) {
  1060.                 continue;
  1061.             }
  1063.             ENa_strand dst_strand = have_strands ?
  1064.                 pseg.GetStrands()[seg*dim + to_row] : eNa_strand_unknown;
  1065.             ENa_strand src_strand = have_strands ?
  1066.                 pseg.GetStrands()[seg*dim + row] : eNa_strand_unknown;
  1068.             // In alignments with multiple sequence types segment length
  1069.             // should be multiplied by character width.
  1070.             TSeqPos src_len = pseg.GetLens()[seg]*src_width*dst_width_rel;
  1071.             TSeqPos dst_len = pseg.GetLens()[seg]*m_Dst_width*src_width_rel;
  1072.             TSeqPos src_start = pseg.GetStarts()[seg*dim + row]*dst_width_rel;
  1073.             TSeqPos dst_start = pseg.GetStarts()[seg*dim + to_row]*src_width_rel;
  1074.             x_NextMappingRange(
  1075.                 src_id, src_start, src_len, src_strand,
  1076.                 dst_id, dst_start, dst_len, dst_strand);
  1077.             _ASSERT(!src_len  &&  !dst_len);
  1078.         }
  1079.     }
  1080. }
  1083. void CSeq_loc_Mapper::x_Initialize(const CSeqMap& seq_map,
  1084.                                    const CSeq_id* top_id)
  1085. {
  1086.     CSeqMap::const_iterator seg_it =
  1087.         seq_map.begin_resolved(m_Scope.GetPointerOrNull(), size_t(-1),
  1088.         CSeqMap::fFindRef);
  1090.     TSeqPos top_start = kInvalidSeqPos;
  1091.     TSeqPos top_stop = kInvalidSeqPos;
  1092.     TSeqPos dst_seg_start = kInvalidSeqPos;
  1093.     CConstRef<CSeq_id> dst_id;
  1095.     for ( ; seg_it; ++seg_it) {
  1096.         _ASSERT(seg_it.GetType() == CSeqMap::eSeqRef);
  1097.         if (seg_it.GetPosition() > top_stop  ||  !dst_id) {
  1098.             // New top-level segment
  1099.             top_start = seg_it.GetPosition();
  1100.             top_stop = seg_it.GetEndPosition() - 1;
  1101.             if (top_id) {
  1102.                 // Top level is a bioseq
  1103.                 dst_id.Reset(top_id);
  1104.                 dst_seg_start = top_start;
  1105.             }
  1106.             else {
  1107.                 // Top level is a seq-loc, positions are
  1108.                 // on the first-level references
  1109.                 dst_id = seg_it.GetRefSeqid().GetSeqId();
  1110.                 dst_seg_start = seg_it.GetRefPosition();
  1111.                 continue;
  1112.             }
  1113.         }
  1114.         // when top_id is set, destination position = GetPosition(),
  1115.         // else it needs to be calculated from top_start/stop and dst_start/stop.
  1116.         TSeqPos dst_from = dst_seg_start + seg_it.GetPosition() - top_start;
  1117.         _ASSERT(dst_from >= dst_seg_start);
  1118.         TSeqPos dst_len = seg_it.GetLength();
  1119.         CConstRef<CSeq_id> src_id(seg_it.GetRefSeqid().GetSeqId());
  1120.         TSeqPos src_from = seg_it.GetRefPosition();
  1121.         TSeqPos src_len = dst_len;
  1122.         ENa_strand src_strand = seg_it.GetRefMinusStrand() ?
  1123.             eNa_strand_minus : eNa_strand_unknown;
  1124.         x_NextMappingRange(*src_id, src_from, src_len, src_strand,
  1125.             *dst_id, dst_from, dst_len, eNa_strand_unknown);
  1126.         _ASSERT(src_len == 0  &&  dst_len == 0);
  1127.     }
  1128. }
  1131. void CSeq_loc_Mapper::x_Initialize(const CSeqMap& seq_map,
  1132.                                    size_t         depth,
  1133.                                    const CSeq_id* top_id)
  1134. {
  1135.     CSeqMap::const_iterator seg_it =
  1136.         seq_map.begin_resolved(m_Scope.GetPointerOrNull(), depth,
  1137.         CSeqMap::fFindRef);
  1139.     TSeqPos top_start = kInvalidSeqPos;
  1140.     TSeqPos top_stop = kInvalidSeqPos;
  1141.     TSeqPos src_seg_start = kInvalidSeqPos;
  1142.     TSeqPos last_stop = kInvalidSeqPos;
  1143.     CConstRef<CSeq_id> src_id;
  1145.     TSeqPos src_from = 0;
  1146.     TSeqPos src_len = 0;
  1147.     TSeqPos dst_from = 0;
  1148.     TSeqPos dst_len = 0;
  1149.     CConstRef<CSeq_id> dst_id;
  1150.     ENa_strand dst_strand = eNa_strand_unknown;
  1152.     for ( ; seg_it; ++seg_it) {
  1153.         _ASSERT(seg_it.GetType() == CSeqMap::eSeqRef);
  1154.         if (seg_it.GetPosition() > top_stop  ||  !src_id) {
  1155.             // New top-level segment
  1156.             top_start = seg_it.GetPosition();
  1157.             top_stop = seg_it.GetEndPosition() - 1;
  1158.             if (top_id) {
  1159.                 // Top level is a bioseq
  1160.                 src_id.Reset(top_id);
  1161.                 src_seg_start = top_start;
  1162.             }
  1163.             else {
  1164.                 // Top level is a seq-loc, positions are
  1165.                 // on the first-level references
  1166.                 src_id = seg_it.GetRefSeqid().GetSeqId();
  1167.                 src_seg_start = seg_it.GetRefPosition();
  1168.                 continue;
  1169.             }
  1170.         }
  1171.         if (seg_it.GetPosition() >= last_stop) {
  1172.             x_NextMappingRange(*src_id, src_from, src_len, eNa_strand_unknown,
  1173.                 *dst_id, dst_from, dst_len, dst_strand);
  1174.             _ASSERT(src_len == 0  &&  dst_len == 0);
  1175.         }
  1176.         src_from = src_seg_start + seg_it.GetPosition() - top_start;
  1177.         src_len = seg_it.GetLength();
  1178.         dst_from = seg_it.GetRefPosition();
  1179.         dst_len = src_len;
  1180.         dst_id.Reset(seg_it.GetRefSeqid().GetSeqId());
  1181.         dst_strand = seg_it.GetRefMinusStrand() ?
  1182.             eNa_strand_minus : eNa_strand_unknown;
  1183.         last_stop = seg_it.GetEndPosition();
  1184.     }
  1185.     if (src_len > 0) {
  1186.         x_NextMappingRange(*src_id, src_from, src_len, eNa_strand_unknown,
  1187.             *dst_id, dst_from, dst_len, dst_strand);
  1188.     }
  1189. }
  1192. CSeq_loc_Mapper::TRangeIterator
  1193. CSeq_loc_Mapper::x_BeginMappingRanges(CSeq_id_Handle id,
  1194.                                       TSeqPos from,
  1195.                                       TSeqPos to)
  1196. {
  1197.     TIdMap::iterator ranges = m_IdMap.find(id);
  1198.     if (ranges == m_IdMap.end()) {
  1199.         return TRangeIterator();
  1200.     }
  1201.     return ranges->second.begin(TRange(from, to));
  1202. }
  1205. bool CSeq_loc_Mapper::x_MapInterval(const CSeq_id&   src_id,
  1206.                                     TRange           src_rg,
  1207.                                     bool             is_set_strand,
  1208.                                     ENa_strand       src_strand,
  1209.                                     TRangeFuzz       orig_fuzz)
  1210. {
  1211.     bool res = false;
  1212.     if (m_UseWidth  &&
  1213.         m_Widths[CSeq_id_Handle::GetHandle(src_id)] & fWidthProtToNuc) {
  1214.         src_rg = TRange(src_rg.GetFrom()*3, src_rg.GetTo()*3 + 2);
  1215.     }
  1217.     CSeq_id_Handle src_idh = CSeq_id_Handle::GetHandle(src_id);
  1219.     // Sorted mapping ranges
  1220.     typedef vector< CRef<CMappingRange> > TSortedMappings;
  1221.     TSortedMappings mappings;
  1222.     TRangeIterator rg_it = x_BeginMappingRanges(
  1223.         src_idh, src_rg.GetFrom(), src_rg.GetTo());
  1224.     for ( ; rg_it; ++rg_it) {
  1225.         mappings.push_back(rg_it->second);
  1226.     }
  1227.     sort(mappings.begin(), mappings.end(), CMappingRangeRef_Less());
  1228.     TSeqPos last_src_to = 0;
  1229.     for (size_t idx = 0; idx < mappings.size(); ++idx) {
  1230.         const CMappingRange& cvt = *mappings[idx];
  1231.         if ( !cvt.CanMap(src_rg.GetFrom(), src_rg.GetTo(),
  1232.             is_set_strand, src_strand) ) {
  1233.             continue;
  1234.         }
  1235.         TSeqPos from = src_rg.GetFrom();
  1236.         TSeqPos to = src_rg.GetTo();
  1237.         TRangeFuzz fuzz = cvt.Map_Fuzz(orig_fuzz);
  1238.         if (from < cvt.m_Src_from) {
  1239.             from = cvt.m_Src_from;
  1240.             // Set partial only if a non-empty range is to be skipped
  1241.             if (cvt.m_Src_from > last_src_to + 1) {
  1242.                 if ( cvt.m_Reverse ) {
  1243.                     if (!fuzz.second) {
  1244.                         fuzz.second.Reset(new CInt_fuzz);
  1245.                     }
  1246.                     fuzz.second->SetLim(CInt_fuzz::eLim_gt);
  1247.                 }
  1248.                 else {
  1249.                     if (!fuzz.first) {
  1250.                         fuzz.first.Reset(new CInt_fuzz);
  1251.                     }
  1252.                     fuzz.first->SetLim(CInt_fuzz::eLim_lt);
  1253.                 }
  1254.                 m_Partial = true;
  1255.             }
  1256.         }
  1257.         last_src_to = cvt.m_Src_to;
  1258.         if ( to > cvt.m_Src_to ) {
  1259.             to = cvt.m_Src_to;
  1260.             // Set partial only if a non-empty range is to be skipped
  1261.             if (idx == mappings.size() - 1 ||
  1262.                 mappings[idx+1]->m_Src_from > last_src_to + 1) {
  1263.                 if ( cvt.m_Reverse ) {
  1264.                     if (!fuzz.first) {
  1265.                         fuzz.first.Reset(new CInt_fuzz);
  1266.                     }
  1267.                     fuzz.first->SetLim(CInt_fuzz::eLim_lt);
  1268.                 }
  1269.                 else {
  1270.                     if (!fuzz.second) {
  1271.                         fuzz.second.Reset(new CInt_fuzz);
  1272.                     }
  1273.                     fuzz.second->SetLim(CInt_fuzz::eLim_gt);
  1274.                 }
  1275.                 m_Partial = true;
  1276.             }
  1277.         }
  1278.         if ( from > to ) {
  1279.             continue;
  1280.         }
  1281.         TRange rg = cvt.Map_Range(from, to);
  1282.         ENa_strand dst_strand;
  1283.         bool is_set_dst_strand = cvt.Map_Strand(is_set_strand,
  1284.             src_strand, &dst_strand);
  1285.         x_GetMappedRanges(cvt.m_Dst_id_Handle,
  1286.             is_set_dst_strand ? int(dst_strand) + 1 : 0)
  1287.             .push_back(TRangeWithFuzz(rg, fuzz));
  1288.         res = true;
  1289.     }
  1290.     return res;
  1291. }
  1294. void CSeq_loc_Mapper::x_PushLocToDstMix(CRef<CSeq_loc> loc)
  1295. {
  1296.     _ASSERT(loc);
  1297.     if ( !m_Dst_loc  ||  !m_Dst_loc->IsMix() ) {
  1298.         CRef<CSeq_loc> tmp = m_Dst_loc;
  1299.         m_Dst_loc.Reset(new CSeq_loc);
  1300.         m_Dst_loc->SetMix();
  1301.         if ( tmp ) {
  1302.             m_Dst_loc->SetMix().Set().push_back(tmp);
  1303.         }
  1304.     }
  1305.     CSeq_loc_mix::Tdata& mix = m_Dst_loc->SetMix().Set();
  1306.     if ( loc->IsNull() ) {
  1307.         if ( m_GapFlag == eGapRemove ) {
  1308.             return;
  1309.         }
  1310.         if ( mix.size() > 0  &&  (*mix.rbegin())->IsNull() ) {
  1311.             // do not create duplicate NULLs
  1312.             return;
  1313.         }
  1314.     }
  1315.     mix.push_back(loc);
  1316. }
  1319. CRef<CSeq_loc> CSeq_loc_Mapper::Map(const CSeq_loc& src_loc)
  1320. {
  1321.     m_Dst_loc.Reset();
  1322.     m_Partial = false; // reset for each location
  1323.     x_MapSeq_loc(src_loc);
  1324.     x_PushRangesToDstMix();
  1325.     x_OptimizeSeq_loc(m_Dst_loc);
  1326.     return m_Dst_loc;
  1327. }
  1330. CRef<CSeq_align> CSeq_loc_Mapper::Map(const CSeq_align& src_align)
  1331. {
  1332.     m_Dst_loc.Reset();
  1333.     m_Partial = false; // reset for each location
  1334.     return x_MapSeq_align(src_align);
  1335. }
  1338. void CSeq_loc_Mapper::x_MapSeq_loc(const CSeq_loc& src_loc)
  1339. {
  1340.     switch ( src_loc.Which() ) {
  1341.     case CSeq_loc::e_Null:
  1342.         if (m_GapFlag == eGapRemove) {
  1343.             return;
  1344.         }
  1345.         // Proceed to seq-loc duplication
  1346.     case CSeq_loc::e_not_set:
  1347.     case CSeq_loc::e_Feat:
  1348.     {
  1349.         x_PushRangesToDstMix();
  1350.         CRef<CSeq_loc> loc(new CSeq_loc);
  1351.         loc->Assign(src_loc);
  1352.         x_PushLocToDstMix(loc);
  1353.         break;
  1354.     }
  1355.     case CSeq_loc::e_Empty:
  1356.     {
  1357.         bool res = false;
  1358.         TRangeIterator mit = x_BeginMappingRanges(
  1359.             CSeq_id_Handle::GetHandle(src_loc.GetEmpty()),
  1360.             TRange::GetWhole().GetFrom(),
  1361.             TRange::GetWhole().GetTo());
  1362.         for ( ; mit; ++mit) {
  1363.             CMappingRange& cvt = *mit->second;
  1364.             if ( cvt.GoodSrcId(src_loc.GetEmpty()) ) {
  1365.                 x_GetMappedRanges(
  1366.                     CSeq_id_Handle::GetHandle(src_loc.GetEmpty()), 0)
  1367.                     .push_back(TRangeWithFuzz(TRange::GetEmpty(),
  1368.                     TRangeFuzz(kEmptyFuzz, kEmptyFuzz)));
  1369.                 res = true;
  1370.                 break;
  1371.             }
  1372.         }
  1373.         if ( !res ) {
  1374.             if ( m_KeepNonmapping ) {
  1375.                 x_PushRangesToDstMix();
  1376.                 CRef<CSeq_loc> loc(new CSeq_loc);
  1377.                 loc->Assign(src_loc);
  1378.                 x_PushLocToDstMix(loc);
  1379.             }
  1380.             else {
  1381.                 m_Partial = true;
  1382.             }
  1383.         }
  1384.         break;
  1385.     }
  1386.     case CSeq_loc::e_Whole:
  1387.     {
  1388.         const CSeq_id& src_id = src_loc.GetWhole();
  1389.         // Convert to the allowed master seq interval
  1390.         TSeqPos src_to = kInvalidSeqPos;
  1391.         if ( m_Scope ) {
  1392.             CBioseq_Handle bh = m_Scope->GetBioseqHandle(src_id);
  1393.             src_to = bh ? bh.GetBioseqLength() : kInvalidSeqPos;
  1394.         }
  1395.         bool res = x_MapInterval(src_id, TRange(0, src_to),
  1396.             false, eNa_strand_unknown,
  1397.             TRangeFuzz(kEmptyFuzz, kEmptyFuzz));
  1398.         if ( !res ) {
  1399.             if ( m_KeepNonmapping ) {
  1400.                 x_PushRangesToDstMix();
  1401.                 CRef<CSeq_loc> loc(new CSeq_loc);
  1402.                 loc->Assign(src_loc);
  1403.                 x_PushLocToDstMix(loc);
  1404.             }
  1405.             else {
  1406.                 m_Partial = true;
  1407.             }
  1408.         }
  1409.         break;
  1410.     }
  1411.     case CSeq_loc::e_Int:
  1412.     {
  1413.         const CSeq_interval& src_int = src_loc.GetInt();
  1414.         TRangeFuzz fuzz(kEmptyFuzz, kEmptyFuzz);
  1415.         if ( src_int.IsSetFuzz_from() ) {
  1416.             fuzz.first.Reset(new CInt_fuzz);
  1417.             fuzz.first->Assign(src_int.GetFuzz_from());
  1418.         }
  1419.         if ( src_int.IsSetFuzz_to() ) {
  1420.             fuzz.second.Reset(new CInt_fuzz);
  1421.             fuzz.second->Assign(src_int.GetFuzz_to());
  1422.         }
  1423.         bool res = x_MapInterval(src_int.GetId(),
  1424.             TRange(src_int.GetFrom(), src_int.GetTo()),
  1425.             src_int.IsSetStrand(),
  1426.             src_int.IsSetStrand() ? src_int.GetStrand() : eNa_strand_unknown,
  1427.             fuzz);
  1428.         if ( !res ) {
  1429.             if ( m_KeepNonmapping ) {
  1430.                 x_PushRangesToDstMix();
  1431.                 CRef<CSeq_loc> loc(new CSeq_loc);
  1432.                 loc->Assign(src_loc);
  1433.                 x_PushLocToDstMix(loc);
  1434.             }
  1435.             else {
  1436.                 m_Partial = true;
  1437.             }
  1438.         }
  1439.         break;
  1440.     }
  1441.     case CSeq_loc::e_Pnt:
  1442.     {
  1443.         const CSeq_point& pnt = src_loc.GetPnt();
  1444.         TRangeFuzz fuzz(kEmptyFuzz, kEmptyFuzz);
  1445.         if ( pnt.IsSetFuzz() ) {
  1446.             fuzz.first.Reset(new CInt_fuzz);
  1447.             fuzz.first->Assign(pnt.GetFuzz());
  1448.         }
  1449.         bool res = x_MapInterval(pnt.GetId(),
  1450.             TRange(pnt.GetPoint(), pnt.GetPoint()),
  1451.             pnt.IsSetStrand(),
  1452.             pnt.IsSetStrand() ? pnt.GetStrand() : eNa_strand_unknown,
  1453.             fuzz);
  1454.         if ( !res ) {
  1455.             if ( m_KeepNonmapping ) {
  1456.                 x_PushRangesToDstMix();
  1457.                 CRef<CSeq_loc> loc(new CSeq_loc);
  1458.                 loc->Assign(src_loc);
  1459.                 x_PushLocToDstMix(loc);
  1460.             }
  1461.             else {
  1462.                 m_Partial = true;
  1463.             }
  1464.         }
  1465.         break;
  1466.     }
  1467.     case CSeq_loc::e_Packed_int:
  1468.     {
  1469.         const CPacked_seqint::Tdata& src_ints = src_loc.GetPacked_int().Get();
  1470.         ITERATE ( CPacked_seqint::Tdata, i, src_ints ) {
  1471.             const CSeq_interval& si = **i;
  1472.             TRangeFuzz fuzz(kEmptyFuzz, kEmptyFuzz);
  1473.             if ( si.IsSetFuzz_from() ) {
  1474.                 fuzz.first.Reset(new CInt_fuzz);
  1475.                 fuzz.first->Assign(si.GetFuzz_from());
  1476.             }
  1477.             if ( si.IsSetFuzz_to() ) {
  1478.                 fuzz.second.Reset(new CInt_fuzz);
  1479.                 fuzz.second->Assign(si.GetFuzz_to());
  1480.             }
  1481.             bool res = x_MapInterval(si.GetId(),
  1482.                 TRange(si.GetFrom(), si.GetTo()),
  1483.                 si.IsSetStrand(),
  1484.                 si.IsSetStrand() ? si.GetStrand() : eNa_strand_unknown,
  1485.                 fuzz);
  1486.             if ( !res ) {
  1487.                 if ( m_KeepNonmapping ) {
  1488.                     x_PushRangesToDstMix();
  1489.                     x_GetMappedRanges(CSeq_id_Handle::GetHandle(si.GetId()),
  1490.                         si.IsSetStrand() ? int(si.GetStrand()) + 1 : 0)
  1491.                         .push_back(TRangeWithFuzz(
  1492.                         TRange(si.GetFrom(), si.GetTo()), fuzz));
  1493.                 }
  1494.                 else {
  1495.                     m_Partial = true;
  1496.                 }
  1497.             }
  1498.         }
  1499.         break;
  1500.     }
  1501.     case CSeq_loc::e_Packed_pnt:
  1502.     {
  1503.         const CPacked_seqpnt& src_pack_pnts = src_loc.GetPacked_pnt();
  1504.         const CPacked_seqpnt::TPoints& src_pnts = src_pack_pnts.GetPoints();
  1505.         ITERATE ( CPacked_seqpnt::TPoints, i, src_pnts ) {
  1506.             TRangeFuzz fuzz(kEmptyFuzz, kEmptyFuzz);
  1507.             if ( src_pack_pnts.IsSetFuzz() ) {
  1508.                 fuzz.first.Reset(new CInt_fuzz);
  1509.                 fuzz.first->Assign(src_pack_pnts.GetFuzz());
  1510.             }
  1511.             bool res = x_MapInterval(
  1512.                 src_pack_pnts.GetId(),
  1513.                 TRange(*i, *i), src_pack_pnts.IsSetStrand(),
  1514.                 src_pack_pnts.IsSetStrand() ?
  1515.                 src_pack_pnts.GetStrand() : eNa_strand_unknown,
  1516.                 fuzz);
  1517.             if ( !res ) {
  1518.                 if ( m_KeepNonmapping ) {
  1519.                     x_PushRangesToDstMix();
  1520.                     x_GetMappedRanges(
  1521.                         CSeq_id_Handle::GetHandle(src_pack_pnts.GetId()),
  1522.                         src_pack_pnts.IsSetStrand() ?
  1523.                         int(src_pack_pnts.GetStrand()) + 1 : 0)
  1524.                         .push_back(TRangeWithFuzz(TRange(*i, *i), fuzz));
  1525.                 }
  1526.                 else {
  1527.                     m_Partial = true;
  1528.                 }
  1529.             }
  1530.         }
  1531.         break;
  1532.     }
  1533.     case CSeq_loc::e_Mix:
  1534.     {
  1535.         x_PushRangesToDstMix();
  1536.         CRef<CSeq_loc> prev = m_Dst_loc;
  1537.         m_Dst_loc.Reset();
  1538.         const CSeq_loc_mix::Tdata& src_mix = src_loc.GetMix().Get();
  1539.         ITERATE ( CSeq_loc_mix::Tdata, i, src_mix ) {
  1540.             x_MapSeq_loc(**i);
  1541.         }
  1542.         x_PushRangesToDstMix();
  1543.         CRef<CSeq_loc> mix = m_Dst_loc;
  1544.         m_Dst_loc = prev;
  1545.         x_OptimizeSeq_loc(mix);
  1546.         x_PushLocToDstMix(mix);
  1547.         break;
  1548.     }
  1549.     case CSeq_loc::e_Equiv:
  1550.     {
  1551.         x_PushRangesToDstMix();
  1552.         CRef<CSeq_loc> prev = m_Dst_loc;
  1553.         m_Dst_loc.Reset();
  1554.         const CSeq_loc_equiv::Tdata& src_equiv = src_loc.GetEquiv().Get();
  1555.         CRef<CSeq_loc> equiv(new CSeq_loc);
  1556.         equiv->SetEquiv();
  1557.         ITERATE ( CSeq_loc_equiv::Tdata, i, src_equiv ) {
  1558.             x_MapSeq_loc(**i);
  1559.             x_PushRangesToDstMix();
  1560.             x_OptimizeSeq_loc(m_Dst_loc);
  1561.             equiv->SetEquiv().Set().push_back(m_Dst_loc);
  1562.             m_Dst_loc.Reset();
  1563.         }
  1564.         m_Dst_loc = prev;
  1565.         x_PushLocToDstMix(equiv);
  1566.         break;
  1567.     }
  1568.     case CSeq_loc::e_Bond:
  1569.     {
  1570.         x_PushRangesToDstMix();
  1571.         CRef<CSeq_loc> prev = m_Dst_loc;
  1572.         m_Dst_loc.Reset();
  1573.         const CSeq_bond& src_bond = src_loc.GetBond();
  1574.         CRef<CSeq_loc> dst_loc(new CSeq_loc);
  1575.         dst_loc->SetBond();
  1576.         CRef<CSeq_loc> pntA;
  1577.         CRef<CSeq_loc> pntB;
  1578.         TRangeFuzz fuzzA(kEmptyFuzz, kEmptyFuzz);
  1579.         if ( src_bond.GetA().IsSetFuzz() ) {
  1580.             fuzzA.first.Reset(new CInt_fuzz);
  1581.             fuzzA.first->Assign(src_bond.GetA().GetFuzz());
  1582.         }
  1583.         bool resA = x_MapInterval(src_bond.GetA().GetId(),
  1584.             TRange(src_bond.GetA().GetPoint(), src_bond.GetA().GetPoint()),
  1585.             src_bond.GetA().IsSetStrand(),
  1586.             src_bond.GetA().IsSetStrand() ?
  1587.             src_bond.GetA().GetStrand() : eNa_strand_unknown,
  1588.             fuzzA);
  1589.         if ( resA ) {
  1590.             pntA = x_GetMappedSeq_loc();
  1591.             _ASSERT(pntA);
  1592.             if ( !pntA->IsPnt() ) {
  1593.                 NCBI_THROW(CLocMapperException, eBadLocation,
  1594.                            "Bond point A can not be mapped to a point");
  1595.             }
  1596.         }
  1597.         else {
  1598.             pntA.Reset(new CSeq_loc);
  1599.             pntA->SetPnt().Assign(src_bond.GetA());
  1600.         }
  1601.         bool resB = false;
  1602.         if ( src_bond.IsSetB() ) {
  1603.             TRangeFuzz fuzzB(kEmptyFuzz, kEmptyFuzz);
  1604.             if ( src_bond.GetB().IsSetFuzz() ) {
  1605.                 fuzzB.first.Reset(new CInt_fuzz);
  1606.                 fuzzB.first->Assign(src_bond.GetB().GetFuzz());
  1607.             }
  1608.             resB = x_MapInterval(src_bond.GetB().GetId(),
  1609.                 TRange(src_bond.GetB().GetPoint(), src_bond.GetB().GetPoint()),
  1610.                 src_bond.GetB().IsSetStrand(),
  1611.                 src_bond.GetB().IsSetStrand() ?
  1612.                 src_bond.GetB().GetStrand() : eNa_strand_unknown,
  1613.                 fuzzB);
  1614.         }
  1615.         if ( resB ) {
  1616.             pntB = x_GetMappedSeq_loc();
  1617.             _ASSERT(pntB);
  1618.             if ( !pntB->IsPnt() ) {
  1619.                 NCBI_THROW(CLocMapperException, eBadLocation,
  1620.                            "Bond point B can not be mapped to a point");
  1621.             }
  1622.         }
  1623.         else {
  1624.             pntB.Reset(new CSeq_loc);
  1625.             pntB->SetPnt().Assign(src_bond.GetB());
  1626.         }
  1627.         m_Dst_loc = prev;
  1628.         if ( resA  ||  resB  ||  m_KeepNonmapping ) {
  1629.             CSeq_bond& dst_bond = dst_loc->SetBond();
  1630.             dst_bond.SetA(pntA->SetPnt());
  1631.             if ( src_bond.IsSetB() ) {
  1632.                 dst_bond.SetB(pntB->SetPnt());
  1633.             }
  1634.             x_PushLocToDstMix(dst_loc);
  1635.         }
  1636.         m_Partial |= (!resA || !resB);
  1637.         break;
  1638.     }
  1639.     default:
  1640.         NCBI_THROW(CAnnotException, eBadLocation,
  1641.                    "Unsupported location type");
  1642.     }
  1643. }
  1646. CRef<CSeq_loc> CSeq_loc_Mapper::x_RangeToSeq_loc(const CSeq_id_Handle& idh,
  1647.                                                  TSeqPos from,
  1648.                                                  TSeqPos to,
  1649.                                                  int strand_idx,
  1650.                                                  TRangeFuzz rg_fuzz)
  1651. {
  1652.     if (m_UseWidth  &&  (m_Widths[idh] & fWidthNucToProt)) {
  1653.         from = from/3;
  1654.         to = to/3;
  1655.     }
  1657.     CRef<CSeq_loc> loc(new CSeq_loc);
  1658.     // If both fuzzes are set, create interval, not point.
  1659.     if (from == to  &&  (!rg_fuzz.first  ||  !rg_fuzz.second)) {
  1660.         // point
  1661.         loc->SetPnt().SetId().Assign(*idh.GetSeqId());
  1662.         loc->SetPnt().SetPoint(from);
  1663.         if (strand_idx > 0) {
  1664.             loc->SetPnt().SetStrand(ENa_strand(strand_idx - 1));
  1665.         }
  1666.         if ( rg_fuzz.first ) {
  1667.             loc->SetPnt().SetFuzz(*rg_fuzz.first);
  1668.         }
  1669.         else if ( rg_fuzz.second ) {
  1670.             loc->SetPnt().SetFuzz(*rg_fuzz.second);
  1671.         }
  1672.     }
  1673.     else {
  1674.         // interval
  1675.         loc->SetInt().SetId().Assign(*idh.GetSeqId());
  1676.         loc->SetInt().SetFrom(from);
  1677.         loc->SetInt().SetTo(to);
  1678.         if (strand_idx > 0) {
  1679.             loc->SetInt().SetStrand(ENa_strand(strand_idx - 1));
  1680.         }
  1681.         if ( rg_fuzz.first ) {
  1682.             loc->SetInt().SetFuzz_from(*rg_fuzz.first);
  1683.         }
  1684.         if ( rg_fuzz.second ) {
  1685.             loc->SetInt().SetFuzz_to(*rg_fuzz.second);
  1686.         }
  1687.     }
  1688.     return loc;
  1689. }
  1692. void CSeq_loc_Mapper::x_OptimizeSeq_loc(CRef<CSeq_loc>& loc)
  1693. {
  1694.     if ( !loc ) {
  1695.         loc.Reset(new CSeq_loc);
  1696.         loc->SetNull();
  1697.         return;
  1698.     }
  1699.     switch (loc->Which()) {
  1700.     case CSeq_loc::e_not_set:
  1701.     case CSeq_loc::e_Feat:
  1702.     case CSeq_loc::e_Null:
  1703.     case CSeq_loc::e_Empty:
  1704.     case CSeq_loc::e_Whole:
  1705.     case CSeq_loc::e_Int:
  1706.     case CSeq_loc::e_Pnt:
  1707.     case CSeq_loc::e_Equiv:
  1708.     case CSeq_loc::e_Bond:
  1709.     case CSeq_loc::e_Packed_int:  // pack ?
  1710.     case CSeq_loc::e_Packed_pnt:  // pack ?
  1711.         return;
  1712.     case CSeq_loc::e_Mix:
  1713.     {
  1714.         switch ( loc->GetMix().Get().size() ) {
  1715.         case 0:
  1716.             loc->SetNull();
  1717.             break;
  1718.         case 1:
  1719.         {
  1720.             CRef<CSeq_loc> single = *loc->SetMix().Set().begin();
  1721.             loc = single;
  1722.             break;
  1723.         }
  1724.         default:
  1725.             break;
  1726.         }
  1727.         break;
  1728.     }
  1729.     default:
  1730.         NCBI_THROW(CAnnotException, eBadLocation,
  1731.                    "Unsupported location type");
  1732.     }
  1733. }
  1736. CRef<CSeq_loc> CSeq_loc_Mapper::x_GetMappedSeq_loc(void)
  1737. {
  1738.     CRef<CSeq_loc> dst_loc(new CSeq_loc);
  1739.     CSeq_loc_mix::Tdata& dst_mix = dst_loc->SetMix().Set();
  1740.     NON_CONST_ITERATE(TRangesById, id_it, m_MappedLocs) {
  1741.         if ( !id_it->first ) {
  1742.             if (m_GapFlag == eGapPreserve) {
  1743.                 CRef<CSeq_loc> null_loc(new CSeq_loc);
  1744.                 null_loc->SetNull();
  1745.                 dst_mix.push_back(null_loc);
  1746.             }
  1747.             continue;
  1748.         }
  1749.         for (int str = 0; str < (int)id_it->second.size(); ++str) {
  1750.             if (id_it->second[str].size() == 0) {
  1751.                 continue;
  1752.             }
  1753.             TSeqPos from = kInvalidSeqPos;
  1754.             TSeqPos to = kInvalidSeqPos;
  1755.             TRangeFuzz fuzz(kEmptyFuzz, kEmptyFuzz);
  1756.             id_it->second[str].sort();
  1757.             NON_CONST_ITERATE(TMappedRanges, rg_it, id_it->second[str]) {
  1758.                 if ( rg_it->first.Empty() ) {
  1759.                     // Empty seq-loc
  1760.                     CRef<CSeq_loc> loc(new CSeq_loc);
  1761.                     loc->SetEmpty().Assign(*id_it->first.GetSeqId());
  1762.                     dst_mix.push_back(loc);
  1763.                     continue;
  1764.                 }
  1765.                 if (to == kInvalidSeqPos) {
  1766.                     from = rg_it->first.GetFrom();
  1767.                     to = rg_it->first.GetTo();
  1768.                     fuzz = rg_it->second;
  1769.                     continue;
  1770.                 }
  1771.                 if (m_MergeFlag == eMergeAbutting) {
  1772.                     if (rg_it->first.GetFrom() == to + 1) {
  1773.                         to = rg_it->first.GetTo();
  1774.                         fuzz.second = rg_it->second.second;
  1775.                         continue;
  1776.                     }
  1777.                 }
  1778.                 if (m_MergeFlag == eMergeAll) {
  1779.                     if (rg_it->first.GetFrom() <= to + 1) {
  1780.                         to = rg_it->first.GetTo();
  1781.                         fuzz.second = rg_it->second.second;
  1782.                         continue;
  1783.                     }
  1784.                 }
  1785.                 // Add new interval or point
  1786.                 dst_mix.push_back(x_RangeToSeq_loc(id_it->first, from, to,
  1787.                     str, fuzz));
  1788.                 from = rg_it->first.GetFrom();
  1789.                 to = rg_it->first.GetTo();
  1790.                 fuzz = rg_it->second;
  1791.             }
  1792.             // last interval or point
  1793.             dst_mix.push_back(x_RangeToSeq_loc(id_it->first, from, to,
  1794.                 str, fuzz));
  1795.         }
  1796.     }
  1797.     m_MappedLocs.clear();
  1798.     x_OptimizeSeq_loc(dst_loc);
  1799.     return dst_loc;
  1800. }
  1803. CRef<CSeq_align> CSeq_loc_Mapper::x_MapSeq_align(const CSeq_align& src_align)
  1804. {
  1805.     CSeq_align_Mapper aln_mapper(src_align);
  1806.     aln_mapper.Convert(*this);
  1807.     return aln_mapper.GetDstAlign();
  1808.     /*
  1809.     ITERATE(TIdMap, id_it, m_IdMap) {
  1810.         ITERATE(TRangeMap, rg_it, id_it->second) {
  1811.             aln_mapper.Convert(*rg_it->second,
  1812.                 m_UseWidth ? m_Widths[id_it->first] : 0);
  1813.         }
  1814.     }
  1815.     */
  1816. }
  1819. END_SCOPE(objects)
  1820. END_NCBI_SCOPE
  1822. /*
  1823. * ---------------------------------------------------------------------------
  1824. * $Log: seq_loc_mapper.cpp,v $
  1825. * Revision 1000.1  2004/06/01 19:24:12  gouriano
  1826. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
  1827. *
  1828. * Revision 1.20  2004/05/26 14:29:20  grichenk
  1829. * Redesigned CSeq_align_Mapper: preserve non-mapping intervals,
  1830. * fixed strands handling, improved performance.
  1831. *
  1832. * Revision 1.19  2004/05/21 21:42:13  gorelenk
  1833. * Added PCH ncbi_pch.hpp
  1834. *
  1835. * Revision 1.18  2004/05/10 20:22:24  grichenk
  1836. * Fixed more warnings (removed inlines)
  1837. *
  1838. * Revision 1.17  2004/05/10 18:26:37  grichenk
  1839. * Fixed 'not used' warnings
  1840. *
  1841. * Revision 1.16  2004/05/07 13:53:18  grichenk
  1842. * Preserve fuzz from original location.
  1843. * Better detection of partial locations.
  1844. *
  1845. * Revision 1.15  2004/05/05 14:04:22  grichenk
  1846. * Use fuzz to indicate truncated intervals. Added KeepNonmapping flag.
  1847. *
  1848. * Revision 1.14  2004/04/23 15:34:49  grichenk
  1849. * Added PreserveDestinationLocs().
  1850. *
  1851. * Revision 1.13  2004/04/12 14:35:59  grichenk
  1852. * Fixed mapping of alignments between nucleotides and proteins
  1853. *
  1854. * Revision 1.12  2004/04/06 13:56:33  grichenk
  1855. * Added possibility to remove gaps (NULLs) from mapped location
  1856. *
  1857. * Revision 1.11  2004/03/31 20:44:26  grichenk
  1858. * Fixed warnings about unhandled cases in switch.
  1859. *
  1860. * Revision 1.10  2004/03/30 21:21:09  grichenk
  1861. * Reduced number of includes.
  1862. *
  1863. * Revision 1.9  2004/03/30 17:00:00  grichenk
  1864. * Fixed warnings, moved inline functions to hpp.
  1865. *
  1866. * Revision 1.8  2004/03/30 15:42:34  grichenk
  1867. * Moved alignment mapper to separate file, added alignment mapping
  1868. * to CSeq_loc_Mapper.
  1869. *
  1870. * Revision 1.7  2004/03/29 15:13:56  grichenk
  1871. * Added mapping down to segments of a bioseq or a seqmap.
  1872. * Fixed: preserve ranges already on the target bioseq(s).
  1873. *
  1874. * Revision 1.6  2004/03/22 21:10:58  grichenk
  1875. * Added mapping from segments to master sequence or through a seq-map.
  1876. *
  1877. * Revision 1.5  2004/03/19 14:19:08  grichenk
  1878. * Added seq-loc mapping through a seq-align.
  1879. *
  1880. * Revision 1.4  2004/03/16 15:47:28  vasilche
  1881. * Added CBioseq_set_Handle and set of EditHandles
  1882. *
  1883. * Revision 1.3  2004/03/11 18:48:02  shomrat
  1884. * skip if empty
  1885. *
  1886. * Revision 1.2  2004/03/11 04:54:48  grichenk
  1887. * Removed inline
  1888. *
  1889. * Revision 1.1  2004/03/10 16:22:29  grichenk
  1890. * Initial revision
  1891. *
  1892. *
  1893. * ===========================================================================
  1894. */