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

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

生物技术

开发平台：

C/C++

seq_align_mapper.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: seq_align_mapper.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 19:23:43  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.10
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: seq_align_mapper.cpp,v 1000.1 2004/06/01 19:23:43 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. *   Alignment mapper
  39. *
  40. */
  42. #include <ncbi_pch.hpp>
  43. #include <objmgr/impl/seq_align_mapper.hpp>
  44. #include <objmgr/seq_id_mapper.hpp>
  45. #include <objmgr/seq_loc_mapper.hpp>
  46. #include <objmgr/objmgr_exception.hpp>
  47. #include <objects/seqalign/seqalign__.hpp>
  48. #include <algorithm>
  50. BEGIN_NCBI_SCOPE
  51. BEGIN_SCOPE(objects)
  54. SAlignment_Segment::SAlignment_Segment(int len, size_t dim)
  55.     : m_Len(len),
  56.       m_Rows(dim),
  57.       m_HaveStrands(false)
  58. {
  59.     return;
  60. }
  63. SAlignment_Segment::SAlignment_Row&
  64. SAlignment_Segment::GetRow(size_t idx)
  65. {
  66.     _ASSERT(m_Rows.size() > idx);
  67.     return m_Rows[idx];
  68. }
  71. SAlignment_Segment::SAlignment_Row&
  72. SAlignment_Segment::CopyRow(size_t idx, const SAlignment_Row& src_row)
  73. {
  74.     SAlignment_Row& dst_row = GetRow(idx);
  75.     dst_row = src_row;
  76.     return dst_row;
  77. }
  80. SAlignment_Segment::SAlignment_Row& SAlignment_Segment::AddRow(size_t idx,
  81.                                                                const CSeq_id& id,
  82.                                                                int start,
  83.                                                                bool is_set_strand,
  84.                                                                ENa_strand strand,
  85.                                                                int width)
  86. {
  87.     SAlignment_Row& row = GetRow(idx);
  88.     row.m_Id = CSeq_id_Handle::GetHandle(id);
  89.     row.m_Start = start < 0 ? kInvalidSeqPos : start;
  90.     row.m_IsSetStrand = is_set_strand;
  91.     row.m_Strand = strand;
  92.     row.m_Width = width;
  93.     m_HaveStrands |= is_set_strand;
  94.     return row;
  95. }
  98. SAlignment_Segment::SAlignment_Row& SAlignment_Segment::AddRow(size_t idx,
  99.                                                                const CSeq_id_Handle& id,
  100.                                                                int start,
  101.                                                                bool is_set_strand,
  102.                                                                ENa_strand strand,
  103.                                                                int width)
  104. {
  105.     SAlignment_Row& row = GetRow(idx);
  106.     row.m_Id = id;
  107.     row.m_Start = start < 0 ? kInvalidSeqPos : start;
  108.     row.m_IsSetStrand = is_set_strand;
  109.     row.m_Strand = strand;
  110.     row.m_Width = width;
  111.     m_HaveStrands |= is_set_strand;
  112.     return row;
  113. }
  116. void SAlignment_Segment::SetScores(const TScores& scores)
  117. {
  118.     m_Scores = scores;
  119. }
  122. CSeq_align_Mapper::CSeq_align_Mapper(const CSeq_align& align)
  123.     : m_OrigAlign(&align),
  124.       m_DstAlign(0),
  125.       m_HaveStrands(false),
  126.       m_HaveWidths(false),
  127.       m_Dim(0),
  128.       m_AlignFlags(eAlign_Normal)
  129. {
  130.     switch ( align.GetSegs().Which() ) {
  131.     case CSeq_align::C_Segs::e_Dendiag:
  132.         x_Init(align.GetSegs().GetDendiag());
  133.         break;
  134.     case CSeq_align::C_Segs::e_Denseg:
  135.         x_Init(align.GetSegs().GetDenseg());
  136.         break;
  137.     case CSeq_align::C_Segs::e_Std:
  138.         x_Init(align.GetSegs().GetStd());
  139.         break;
  140.     case CSeq_align::C_Segs::e_Packed:
  141.         x_Init(align.GetSegs().GetPacked());
  142.         break;
  143.     case CSeq_align::C_Segs::e_Disc:
  144.         x_Init(align.GetSegs().GetDisc());
  145.         break;
  146.     default:
  147.         break;
  148.     }
  149. }
  152. SAlignment_Segment& CSeq_align_Mapper::x_PushSeg(int len, size_t dim)
  153. {
  154.     m_Segs.push_back(SAlignment_Segment(len, dim));
  155.     return m_Segs.back();
  156. }
  159. SAlignment_Segment& CSeq_align_Mapper::x_InsertSeg(TSegments::iterator& where,
  160.                                                    int len,
  161.                                                    size_t dim)
  162. {
  163.     return *m_Segs.insert(where, SAlignment_Segment(len, dim));
  164. }
  167. void CSeq_align_Mapper::x_Init(const TDendiag& diags)
  168. {
  169.     ITERATE(TDendiag, diag_it, diags) {
  170.         const CDense_diag& diag = **diag_it;
  171.         size_t dim = diag.GetDim();
  172.         if (dim != diag.GetIds().size()) {
  173.             ERR_POST(Warning << "Invalid 'ids' size in dendiag");
  174.             dim = min(dim, diag.GetIds().size());
  175.         }
  176.         if (dim != diag.GetStarts().size()) {
  177.             ERR_POST(Warning << "Invalid 'starts' size in dendiag");
  178.             dim = min(dim, diag.GetStarts().size());
  179.         }
  180.         m_HaveStrands = diag.IsSetStrands();
  181.         if (m_HaveStrands && dim != diag.GetStrands().size()) {
  182.             ERR_POST(Warning << "Invalid 'strands' size in dendiag");
  183.             dim = min(dim, diag.GetStrands().size());
  184.         }
  185.         if (dim != m_Dim) {
  186.             if ( m_Dim ) {
  187.                 m_AlignFlags = eAlign_MultiDim;
  188.             }
  189.             m_Dim = max(dim, m_Dim);
  190.         }
  191.         SAlignment_Segment& seg = x_PushSeg(diag.GetLen(), dim);
  192.         ENa_strand strand = eNa_strand_unknown;
  193.         if ( diag.IsSetScores() ) {
  194.             seg.SetScores(diag.GetScores());
  195.         }
  196.         for (size_t row = 0; row < dim; ++row) {
  197.             if ( m_HaveStrands ) {
  198.                 strand = diag.GetStrands()[row];
  199.             }
  200.             seg.AddRow(row,
  201.                 *diag.GetIds()[row],
  202.                 diag.GetStarts()[row],
  203.                 m_HaveStrands,
  204.                 strand,
  205.                 0);
  206.         }
  207.     }
  208. }
  211. void CSeq_align_Mapper::x_Init(const CDense_seg& denseg)
  212. {
  213.     m_Dim = denseg.GetDim();
  214.     size_t numseg = denseg.GetNumseg();
  215.     // claimed dimension may not be accurate :-/
  216.     if (numseg != denseg.GetLens().size()) {
  217.         ERR_POST(Warning << "Invalid 'lens' size in denseg");
  218.         numseg = min(numseg, denseg.GetLens().size());
  219.     }
  220.     if (m_Dim != denseg.GetIds().size()) {
  221.         ERR_POST(Warning << "Invalid 'ids' size in denseg");
  222.         m_Dim = min(m_Dim, denseg.GetIds().size());
  223.     }
  224.     if (m_Dim*numseg != denseg.GetStarts().size()) {
  225.         ERR_POST(Warning << "Invalid 'starts' size in denseg");
  226.         m_Dim = min(m_Dim*numseg, denseg.GetStarts().size()) / numseg;
  227.     }
  228.     m_HaveStrands = denseg.IsSetStrands();
  229.     if (m_HaveStrands && m_Dim*numseg != denseg.GetStrands().size()) {
  230.         ERR_POST(Warning << "Invalid 'strands' size in denseg");
  231.         m_Dim = min(m_Dim*numseg, denseg.GetStrands().size()) / numseg;
  232.     }
  233.     m_HaveWidths = denseg.IsSetWidths();
  234.     ENa_strand strand = eNa_strand_unknown;
  235.     for (size_t seg = 0;  seg < numseg;  seg++) {
  236.         SAlignment_Segment& alnseg = x_PushSeg(denseg.GetLens()[seg], m_Dim);
  237.         if ( seg == 0  &&  denseg.IsSetScores() ) {
  238.             // Set scores for the first segment only to avoid multiple copies
  239.             alnseg.SetScores(denseg.GetScores());
  240.         }
  241.         for (unsigned int row = 0;  row < m_Dim;  row++) {
  242.             if ( m_HaveStrands ) {
  243.                 strand = denseg.GetStrands()[seg*m_Dim + row];
  244.             }
  245.             alnseg.AddRow(row,
  246.                 *denseg.GetIds()[row],
  247.                 denseg.GetStarts()[seg*m_Dim + row],
  248.                 m_HaveStrands,
  249.                 strand,
  250.                 m_HaveWidths ? denseg.GetWidths()[row] : 0);
  251.         }
  252.     }
  253. }
  256. void CSeq_align_Mapper::x_Init(const TStd& sseg)
  257. {
  258.     typedef map<CSeq_id_Handle, int> TSegLenMap;
  259.     TSegLenMap seglens;
  261.     ITERATE ( CSeq_align::C_Segs::TStd, it, sseg ) {
  262.         const CStd_seg& stdseg = **it;
  263.         size_t dim = stdseg.GetDim();
  264.         if (stdseg.IsSetIds()
  265.             && dim != stdseg.GetIds().size()) {
  266.             ERR_POST(Warning << "Invalid 'ids' size in std-seg");
  267.             dim = min(dim, stdseg.GetIds().size());
  268.         }
  269.         SAlignment_Segment& seg = x_PushSeg(0, dim);
  270.         if ( stdseg.IsSetScores() ) {
  271.             seg.SetScores(stdseg.GetScores());
  272.         }
  273.         int seg_len = 0;
  274.         bool multi_width = false;
  275.         ITERATE ( CStd_seg::TLoc, it_loc, (*it)->GetLoc() ) {
  276.             const CSeq_loc& loc = **it_loc;
  277.             const CSeq_id* id = 0;
  278.             int start = -1;
  279.             int len = 0;
  280.             ENa_strand strand = eNa_strand_unknown;
  281.             bool have_strand = false;
  282.             unsigned int row_idx = 0;
  283.             for (CSeq_loc_CI rg(loc); rg; ++rg, ++row_idx) {
  284.                 if (row_idx > dim) {
  285.                     ERR_POST(Warning << "Invalid number of rows in std-seg");
  286.                     dim = row_idx;
  287.                     seg.m_Rows.resize(dim);
  288.                 }
  289.                 id = &rg.GetSeq_id();
  290.                 if ( rg.IsEmpty() ) {
  291.                     start = kInvalidSeqPos;
  292.                 }
  293.                 else {
  294.                     start = rg.GetRange().GetFrom();
  295.                     len = rg.GetRange().GetLength();
  296.                     strand = rg.GetStrand();
  297.                     if (strand != eNa_strand_unknown) {
  298.                         m_HaveStrands = have_strand = true;
  299.                     }
  300.                 }
  301.                 if (len > 0) {
  302.                     if (seg_len == 0) {
  303.                         seg_len = len;
  304.                     }
  305.                     else if (len != seg_len) {
  306.                         multi_width = true;
  307.                         if (len/3 == seg_len) {
  308.                             seg_len = len;
  309.                         }
  310.                         else if (len*3 != seg_len) {
  311.                             NCBI_THROW(CAnnotException, eBadLocation,
  312.                                     "Rows have different lengths in std-seg");
  313.                         }
  314.                     }
  315.                 }
  316.                 seglens[seg.AddRow(row_idx,
  317.                     *id,
  318.                     start,
  319.                     have_strand,
  320.                     strand,
  321.                     0).m_Id] = len;
  322.             }
  323.             if (dim != m_Dim) {
  324.                 if ( m_Dim ) {
  325.                     m_AlignFlags = eAlign_MultiDim;
  326.                 }
  327.                 m_Dim = max(dim, m_Dim);
  328.             }
  329.         }
  330.         seg.m_Len = seg_len;
  331.         if ( multi_width ) {
  332.             // Adjust each segment width. Do not check if sequence always
  333.             // has the same width.
  334.             for (size_t i = 0; i < seg.m_Rows.size(); ++i) {
  335.                 if (seglens[seg.m_Rows[i].m_Id] != seg_len) {
  336.                     seg.m_Rows[i].m_Width = 3;
  337.                 }
  338.             }
  339.         }
  340.         seglens.clear();
  341.         m_HaveWidths |= multi_width;
  342.     }
  343. }
  346. void CSeq_align_Mapper::x_Init(const CPacked_seg& pseg)
  347. {
  348.     m_Dim = pseg.GetDim();
  349.     size_t numseg = pseg.GetNumseg();
  350.     // claimed dimension may not be accurate :-/
  351.     if (numseg != pseg.GetLens().size()) {
  352.         ERR_POST(Warning << "Invalid 'lens' size in packed-seg");
  353.         numseg = min(numseg, pseg.GetLens().size());
  354.     }
  355.     if (m_Dim != pseg.GetIds().size()) {
  356.         ERR_POST(Warning << "Invalid 'ids' size in packed-seg");
  357.         m_Dim = min(m_Dim, pseg.GetIds().size());
  358.     }
  359.     if (m_Dim*numseg != pseg.GetStarts().size()) {
  360.         ERR_POST(Warning << "Invalid 'starts' size in packed-seg");
  361.         m_Dim = min(m_Dim*numseg, pseg.GetStarts().size()) / numseg;
  362.     }
  363.     m_HaveStrands = pseg.IsSetStrands();
  364.     if (m_HaveStrands && m_Dim*numseg != pseg.GetStrands().size()) {
  365.         ERR_POST(Warning << "Invalid 'strands' size in packed-seg");
  366.         m_Dim = min(m_Dim*numseg, pseg.GetStrands().size()) / numseg;
  367.     }
  368.     ENa_strand strand = eNa_strand_unknown;
  369.     for (size_t seg = 0;  seg < numseg;  seg++) {
  370.         SAlignment_Segment& alnseg = x_PushSeg(pseg.GetLens()[seg], m_Dim);
  371.         if ( seg == 0  &&  pseg.IsSetScores() ) {
  372.             // Set scores for the first segment only to avoid multiple copies
  373.             alnseg.SetScores(pseg.GetScores());
  374.         }
  375.         for (unsigned int row = 0;  row < m_Dim;  row++) {
  376.             if ( m_HaveStrands ) {
  377.                 strand = pseg.GetStrands()[seg*m_Dim + row];
  378.             }
  379.             alnseg.AddRow(row, 
  380.                 *pseg.GetIds()[row],
  381.                 (pseg.GetPresent()[seg*m_Dim + row] ?
  382.                 pseg.GetStarts()[seg*m_Dim + row] : kInvalidSeqPos),
  383.                 m_HaveStrands,
  384.                 strand,
  385.                 0);
  386.         }
  387.     }
  388. }
  391. void CSeq_align_Mapper::x_Init(const CSeq_align_set& align_set)
  392. {
  393.     const CSeq_align_set::Tdata& data = align_set.Get();
  394.     ITERATE(CSeq_align_set::Tdata, it, data) {
  395.         m_SubAligns.push_back(CSeq_align_Mapper(**it));
  396.     }
  397. }
  400. // Mapping through CSeq_loc_Conversion_Set
  402. struct CConversionRef_Less
  403. {
  404.     bool operator()(const CRef<CSeq_loc_Conversion>& x,
  405.                     const CRef<CSeq_loc_Conversion>& y) const;
  406. };
  409. bool CConversionRef_Less::operator()(const CRef<CSeq_loc_Conversion>& x,
  410.                                      const CRef<CSeq_loc_Conversion>& y) const
  411. {
  412.     if (x->m_Src_id_Handle != y->m_Src_id_Handle) {
  413.         return x->m_Src_id_Handle < y->m_Src_id_Handle;
  414.     }
  415.     // Leftmost first
  416.     if (x->m_Src_from != y->m_Src_from) {
  417.         return x->m_Src_from < y->m_Src_from;
  418.     }
  419.     // Longest first
  420.     return x->m_Src_to > y->m_Src_to;
  421. }
  424. void CSeq_align_Mapper::Convert(CSeq_loc_Conversion_Set& cvts,
  425.                                 unsigned int loc_index_shift)
  426. {
  427.     m_DstAlign.Reset();
  429.     if (m_SubAligns.size() > 0) {
  430.         NON_CONST_ITERATE(TSubAligns, it, m_SubAligns) {
  431.             it->Convert(cvts, loc_index_shift);
  432.             loc_index_shift += it->m_Dim;
  433.         }
  434.         return;
  435.     }
  436.     x_ConvertAlign(cvts, loc_index_shift);
  437. }
  440. void CSeq_align_Mapper::x_ConvertAlign(CSeq_loc_Conversion_Set& cvts,
  441.                                        unsigned int loc_index_shift)
  442. {
  443.     if (cvts.m_CvtByIndex.size() == 0) {
  444.         // Single mapping
  445.         _ASSERT(cvts.m_SingleConv);
  446.         x_ConvertRow(*cvts.m_SingleConv,
  447.             cvts.m_SingleIndex - loc_index_shift);
  448.         return;
  449.     }
  450.     CSeq_loc_Conversion_Set::TConvByIndex::iterator idx_it =
  451.         cvts.m_CvtByIndex.lower_bound(loc_index_shift);
  452.     for ( ; idx_it != cvts.m_CvtByIndex.end()
  453.         &&  idx_it->first < loc_index_shift + m_Dim; ++idx_it) {
  454.         x_ConvertRow(idx_it->second, idx_it->first - loc_index_shift);
  455.     }
  456. }
  459. void CSeq_align_Mapper::x_ConvertRow(CSeq_loc_Conversion& cvt,
  460.                                      size_t row)
  461. {
  462.     CSeq_id_Handle dst_id;
  463.     TSegments::iterator seg_it = m_Segs.begin();
  464.     for ( ; seg_it != m_Segs.end(); ) {
  465.         if (seg_it->m_Rows.size() <= row) {
  466.             // No such row in the current segment
  467.             ++seg_it;
  468.             m_AlignFlags = eAlign_MultiDim;
  469.             continue;
  470.         }
  471.         CSeq_id_Handle seg_id = x_ConvertSegment(seg_it, cvt, row);
  472.         if (dst_id) {
  473.             if (dst_id != seg_id  &&  m_AlignFlags == eAlign_Normal) {
  474.                 m_AlignFlags = eAlign_MultiId;
  475.             }
  476.             dst_id = seg_id;
  477.         }
  478.     }
  479. }
  482. void CSeq_align_Mapper::x_ConvertRow(TIdMap& cvts,
  483.                                      size_t row)
  484. {
  485.     CSeq_id_Handle dst_id;
  486.     TSegments::iterator seg_it = m_Segs.begin();
  487.     for ( ; seg_it != m_Segs.end(); ) {
  488.         if (seg_it->m_Rows.size() <= row) {
  489.             // No such row in the current segment
  490.             ++seg_it;
  491.             m_AlignFlags = eAlign_MultiDim;
  492.             continue;
  493.         }
  494.         CSeq_id_Handle seg_id = x_ConvertSegment(seg_it, cvts, row);
  495.         if (dst_id) {
  496.             if (dst_id != seg_id  &&  m_AlignFlags == eAlign_Normal) {
  497.                 m_AlignFlags = eAlign_MultiId;
  498.             }
  499.             dst_id = seg_id;
  500.         }
  501.     }
  502. }
  505. CSeq_id_Handle
  506. CSeq_align_Mapper::x_ConvertSegment(TSegments::iterator& seg_it,
  507.                                     CSeq_loc_Conversion& cvt,
  508.                                     size_t row)
  509. {
  510.     TSegments::iterator old_it = seg_it;
  511.     ++seg_it;
  512.     SAlignment_Segment::SAlignment_Row& aln_row = old_it->m_Rows[row];
  513.     if (aln_row.m_Id != cvt.m_Src_id_Handle) {
  514.         return aln_row.m_Id;
  515.     }
  516.     if (aln_row.m_Start == kInvalidSeqPos) {
  517.         // ??? Skipped row - change ID
  518.         aln_row.m_Id = cvt.m_Dst_id_Handle;
  519.         aln_row.SetMapped();
  520.         return aln_row.m_Id;
  521.     }
  522.     TRange rg(aln_row.m_Start, aln_row.m_Start + old_it->m_Len - 1);
  523.     if (!cvt.ConvertInterval(rg.GetFrom(), rg.GetTo(), aln_row.m_Strand) ) {
  524.         // Do not erase the segment, just change the row ID and reset start
  525.         aln_row.m_Start = kInvalidSeqPos;
  526.         aln_row.m_Id = cvt.m_Dst_id_Handle;
  527.         aln_row.SetMapped();
  528.         return aln_row.m_Id;
  529.     }
  531.     // At least part of the interval was converted.
  532.     TSeqPos dl = cvt.m_Src_from <= rg.GetFrom() ?
  533.         0 : cvt.m_Src_from - rg.GetFrom();
  534.     TSeqPos dr = cvt.m_Src_to >= rg.GetTo() ?
  535.         0 : rg.GetTo() - cvt.m_Src_to;
  536.     if (dl > 0) {
  537.         // Add segment for the skipped range
  538.         SAlignment_Segment& lseg = x_InsertSeg(seg_it, dl,
  539.             old_it->m_Rows.size());
  540.         for (size_t r = 0; r < old_it->m_Rows.size(); ++r) {
  541.             SAlignment_Segment::SAlignment_Row& lrow =
  542.                 lseg.CopyRow(r, old_it->m_Rows[r]);
  543.             if (r == row) {
  544.                 lrow.m_Start = kInvalidSeqPos;
  545.                 lrow.m_Id = cvt.m_Dst_id_Handle;
  546.             }
  547.             else if (lrow.m_Start != kInvalidSeqPos &&
  548.                 !lrow.SameStrand(aln_row)) {
  549.                 // Adjust start for minus strand
  550.                 lrow.m_Start += old_it->m_Len - lseg.m_Len;
  551.             }
  552.         }
  553.     }
  554.     rg.SetFrom(rg.GetFrom() + dl);
  555.     SAlignment_Segment& mseg = x_InsertSeg(seg_it,
  556.         rg.GetLength() - dr,
  557.         old_it->m_Rows.size());
  558.     if (!dl  &&  !dr) {
  559.         // copy scores if there's no truncation
  560.         mseg.SetScores(old_it->m_Scores);
  561.     }
  562.     for (size_t r = 0; r < old_it->m_Rows.size(); ++r) {
  563.         SAlignment_Segment::SAlignment_Row& mrow =
  564.             mseg.CopyRow(r, old_it->m_Rows[r]);
  565.         if (r == row) {
  566.             // translate id and coords
  567.             mrow.m_Id = cvt.m_Dst_id_Handle;
  568.             mrow.m_Start = cvt.m_LastRange.GetFrom();
  569.             mrow.m_IsSetStrand |= (cvt.m_LastStrand != eNa_strand_unknown);
  570.             mrow.m_Strand = cvt.m_LastStrand;
  571.             mrow.SetMapped();
  572.             mseg.m_HaveStrands |= mrow.m_IsSetStrand;
  573.         }
  574.         else {
  575.             if (mrow.m_Start != kInvalidSeqPos) {
  576.                 if (mrow.SameStrand(aln_row)) {
  577.                     mrow.m_Start += dl;
  578.                 }
  579.                 else {
  580.                     mrow.m_Start += old_it->m_Len - dl - mseg.m_Len;
  581.                 }
  582.             }
  583.         }
  584.     }
  585.     cvt.m_LastType = cvt.eMappedObjType_not_set;
  586.     dl += rg.GetLength() - dr;
  587.     rg.SetFrom(rg.GetTo() - dr);
  588.     if (dr > 0) {
  589.         // Add the remaining unmapped range
  590.         SAlignment_Segment& rseg = x_InsertSeg(seg_it,
  591.             dr,
  592.             old_it->m_Rows.size());
  593.         for (size_t r = 0; r < old_it->m_Rows.size(); ++r) {
  594.             SAlignment_Segment::SAlignment_Row& rrow =
  595.                 rseg.CopyRow(r, old_it->m_Rows[r]);
  596.             if (r == row) {
  597.                 rrow.m_Start = kInvalidSeqPos;
  598.                 rrow.m_Id = cvt.m_Dst_id_Handle;
  599.             }
  600.             else {
  601.                 if (rrow.SameStrand(aln_row)) {
  602.                     rrow.m_Start += dl;
  603.                 }
  604.             }
  605.         }
  606.     }
  607.     m_Segs.erase(old_it);
  608.     return cvt.m_Dst_id_Handle;
  609. }
  612. CSeq_id_Handle
  613. CSeq_align_Mapper::x_ConvertSegment(TSegments::iterator& seg_it,
  614.                                     TIdMap& id_map,
  615.                                     size_t row)
  616. {
  617.     TSegments::iterator old_it = seg_it;
  618.     SAlignment_Segment& seg = *old_it;
  619.     ++seg_it;
  620.     SAlignment_Segment::SAlignment_Row& aln_row = seg.m_Rows[row];
  621.     if (aln_row.m_Start == kInvalidSeqPos) {
  622.         // skipped row
  623.         return aln_row.m_Id;
  624.     }
  625.     TRange rg(aln_row.m_Start, aln_row.m_Start + seg.m_Len - 1);
  626.     TIdMap::iterator id_it = id_map.find(aln_row.m_Id);
  627.     if (id_it == id_map.end()) {
  628.         // ID not found in the segment, leave the row unchanged
  629.         return aln_row.m_Id;
  630.     }
  631.     TRangeMap& rmap = id_it->second;
  632.     if ( rmap.empty() ) {
  633.         // No mappings for this segment
  634.         return aln_row.m_Id;
  635.     }
  636.     // Sorted mappings
  637.     typedef vector< CRef<CSeq_loc_Conversion> > TSortedConversions;
  638.     TSortedConversions cvts;
  639.     for (TRangeMap::iterator rg_it = rmap.begin(rg); rg_it; ++rg_it) {
  640.         cvts.push_back(rg_it->second);
  641.     }
  642.     sort(cvts.begin(), cvts.end(), CConversionRef_Less());
  644.     bool mapped = false;
  645.     CSeq_id_Handle dst_id;
  646.     EAlignFlags align_flags = eAlign_Normal;
  647.     TSeqPos left_shift = 0;
  648.     for (size_t cvt_idx = 0; cvt_idx < cvts.size(); ++cvt_idx) {
  649.         CSeq_loc_Conversion& cvt = *cvts[cvt_idx];
  650.         if (!cvt.ConvertInterval(rg.GetFrom(), rg.GetTo(), aln_row.m_Strand) ) {
  651.             continue;
  652.         }
  653.         // Check destination id
  654.         if ( dst_id ) {
  655.             if (cvt.m_Dst_id_Handle != dst_id) {
  656.                 align_flags = eAlign_MultiId;
  657.             }
  658.         }
  659.         dst_id = cvt.m_Dst_id_Handle;
  661.         // At least part of the interval was converted.
  662.         TSeqPos dl = cvt.m_Src_from <= rg.GetFrom() ?
  663.             0 : cvt.m_Src_from - rg.GetFrom();
  664.         TSeqPos dr = cvt.m_Src_to >= rg.GetTo() ?
  665.             0 : rg.GetTo() - cvt.m_Src_to;
  666.         if (dl > 0) {
  667.             // Add segment for the skipped range
  668.             SAlignment_Segment& lseg = x_InsertSeg(seg_it,
  669.                 dl, seg.m_Rows.size());
  670.             for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  671.                 SAlignment_Segment::SAlignment_Row& lrow =
  672.                     lseg.CopyRow(r, seg.m_Rows[r]);
  673.                 if (r == row) {
  674.                     lrow.m_Start = kInvalidSeqPos;
  675.                     lrow.m_Id = dst_id;
  676.                 }
  677.                 else if (lrow.m_Start != kInvalidSeqPos) {
  678.                     if (lrow.SameStrand(aln_row)) {
  679.                         lrow.m_Start += left_shift;
  680.                     }
  681.                     else {
  682.                         lrow.m_Start += seg.m_Len - lseg.m_Len - left_shift;
  683.                     }
  684.                 }
  685.             }
  686.         }
  687.         left_shift += dl;
  688.         SAlignment_Segment& mseg = x_InsertSeg(seg_it,
  689.             rg.GetLength() - dl - dr, seg.m_Rows.size());
  690.         if (!dl  &&  !dr) {
  691.             // copy scores if there's no truncation
  692.             mseg.SetScores(seg.m_Scores);
  693.         }
  694.         for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  695.             SAlignment_Segment::SAlignment_Row& mrow =
  696.                 mseg.CopyRow(r, seg.m_Rows[r]);
  697.             if (r == row) {
  698.                 // translate id and coords
  699.                 mrow.m_Id = cvt.m_Dst_id_Handle;
  700.                 mrow.m_Start = cvt.m_LastRange.GetFrom();
  701.                 mrow.m_IsSetStrand |= (cvt.m_LastStrand != eNa_strand_unknown);
  702.                 mrow.m_Strand = cvt.m_LastStrand;
  703.                 mrow.SetMapped();
  704.                 mseg.m_HaveStrands |= mrow.m_IsSetStrand;
  705.             }
  706.             else {
  707.                 if (mrow.m_Start != kInvalidSeqPos) {
  708.                     if (mrow.SameStrand(aln_row)) {
  709.                         mrow.m_Start += left_shift;
  710.                     }
  711.                     else {
  712.                         mrow.m_Start += seg.m_Len - left_shift - mseg.m_Len;
  713.                     }
  714.                 }
  715.             }
  716.         }
  717.         cvt.m_LastType = cvt.eMappedObjType_not_set;
  718.         mapped = true;
  719.         left_shift += mseg.m_Len;
  720.         rg.SetFrom(aln_row.m_Start + left_shift);
  721.     }
  722.     if (align_flags == eAlign_MultiId  &&  m_AlignFlags == eAlign_Normal) {
  723.         m_AlignFlags = align_flags;
  724.     }
  725.     if ( !mapped ) {
  726.         // Do not erase the segment, just change the row ID and reset start
  727.         seg.m_Rows[row].m_Start = kInvalidSeqPos;
  728.         seg.m_Rows[row].m_Id = rmap.begin()->second->m_Dst_id_Handle;
  729.         seg.m_Rows[row].SetMapped();
  730.         return seg.m_Rows[row].m_Id;
  731.     }
  732.     if (rg.GetFrom() <= rg.GetTo()) {
  733.         // Add the remaining unmapped range
  734.         SAlignment_Segment& rseg = x_InsertSeg(seg_it,
  735.             rg.GetLength(), seg.m_Rows.size());
  736.         for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  737.             SAlignment_Segment::SAlignment_Row& rrow =
  738.                 rseg.CopyRow(r, seg.m_Rows[r]);
  739.             if (r == row) {
  740.                 rrow.m_Start = kInvalidSeqPos;
  741.                 rrow.m_Id = dst_id;
  742.             }
  743.             else if (rrow.m_Start != kInvalidSeqPos) {
  744.                 if (rrow.SameStrand(aln_row)) {
  745.                     rrow.m_Start += left_shift;
  746.                 }
  747.             }
  748.         }
  749.     }
  750.     m_Segs.erase(old_it);
  751.     return align_flags == eAlign_MultiId ? CSeq_id_Handle() : dst_id;
  752. }
  755. // Mapping through CSeq_loc_Mapper
  757. void CSeq_align_Mapper::Convert(CSeq_loc_Mapper& mapper)
  758. {
  759.     m_DstAlign.Reset();
  761.     if (m_SubAligns.size() > 0) {
  762.         NON_CONST_ITERATE(TSubAligns, it, m_SubAligns) {
  763.             it->Convert(mapper);
  764.         }
  765.         return;
  766.     }
  767.     x_ConvertAlign(mapper);
  768. }
  771. void CSeq_align_Mapper::x_ConvertAlign(CSeq_loc_Mapper& mapper)
  772. {
  773.     if (m_Segs.size() == 0) {
  774.         return;
  775.     }
  776.     for (size_t row_idx = 0; row_idx < m_Dim; ++row_idx) {
  777.         x_ConvertRow(mapper, row_idx);
  778.     }
  779. }
  782. void CSeq_align_Mapper::x_ConvertRow(CSeq_loc_Mapper& mapper,
  783.                                      size_t row)
  784. {
  785.     CSeq_id_Handle dst_id;
  786.     TSegments::iterator seg_it = m_Segs.begin();
  787.     for ( ; seg_it != m_Segs.end(); ) {
  788.         if (seg_it->m_Rows.size() <= row) {
  789.             // No such row in the current segment
  790.             ++seg_it;
  791.             m_AlignFlags = eAlign_MultiDim;
  792.             continue;
  793.         }
  794.         CSeq_id_Handle seg_id = x_ConvertSegment(seg_it, mapper, row);
  795.         if (dst_id) {
  796.             if (dst_id != seg_id  &&  m_AlignFlags == eAlign_Normal) {
  797.                 m_AlignFlags = eAlign_MultiId;
  798.             }
  799.             dst_id = seg_id;
  800.         }
  801.     }
  802. }
  805. CSeq_id_Handle
  806. CSeq_align_Mapper::x_ConvertSegment(TSegments::iterator& seg_it,
  807.                                     CSeq_loc_Mapper& mapper,
  808.                                     size_t row)
  809. {
  810.     TSegments::iterator old_it = seg_it;
  811.     SAlignment_Segment& seg = *old_it;
  812.     ++seg_it;
  813.     SAlignment_Segment::SAlignment_Row& aln_row = seg.m_Rows[row];
  814.     if (aln_row.m_Start == kInvalidSeqPos) {
  815.         // skipped row
  816.         return aln_row.m_Id;
  817.     }
  818.     CSeq_loc_Mapper::TIdMap::iterator id_it =
  819.         mapper.m_IdMap.find(aln_row.m_Id);
  820.     if (id_it == mapper.m_IdMap.end()) {
  821.         // ID not found in the segment, leave the row unchanged
  822.         return aln_row.m_Id;
  823.     }
  824.     CSeq_loc_Mapper::TRangeMap& rmap = id_it->second;
  825.     if ( rmap.empty() ) {
  826.         // No mappings for this segment
  827.         return aln_row.m_Id;
  828.     }
  829.     // Sorted mappings
  830.     typedef vector< CRef<CMappingRange> > TSortedMappings;
  831.     TSortedMappings mappings;
  832.     CSeq_loc_Mapper::TRangeMap::iterator rg_it = rmap.begin();
  833.     for ( ; rg_it; ++rg_it) {
  834.         mappings.push_back(rg_it->second);
  835.     }
  836.     sort(mappings.begin(), mappings.end(), CMappingRangeRef_Less());
  838.     bool mapped = false;
  839.     CSeq_id_Handle dst_id;
  840.     EAlignFlags align_flags = eAlign_Normal;
  841.     int width_flag = mapper.m_UseWidth ?
  842.         mapper.m_Widths[aln_row.m_Id] : 0;
  843.     int dst_width = width_flag ? 3 : 1;
  844.     TSeqPos start = aln_row.m_Start;
  845.     if (width_flag & CSeq_loc_Mapper::fWidthProtToNuc) {
  846.         dst_width = 1;
  847.         if (aln_row.m_Width == 3) {
  848.             start *= 3;
  849.         }
  850.     }
  851.     else if (!width_flag  &&  aln_row.m_Width == 3) {
  852.         dst_width = 3;
  853.         start *= 3;
  854.     }
  855.     TSeqPos stop = start + seg.m_Len - 1;
  856.     TSeqPos left_shift = 0;
  857.     for (size_t map_idx = 0; map_idx < mappings.size(); ++map_idx) {
  858.         CRef<CMappingRange> mapping(mappings[map_idx]);
  859.         // Adjust mapping coordinates according to width
  860.         if (width_flag & CSeq_loc_Mapper::fWidthProtToNuc) {
  861.             if (width_flag & CSeq_loc_Mapper::fWidthNucToProt) {
  862.                 // Copy mapping
  863.                 mapping.Reset(new CMappingRange(*mappings[map_idx]));
  864.                 mapping->m_Src_from *= 3;
  865.                 mapping->m_Src_to = (mapping->m_Src_to + 1)*3 - 1;
  866.                 mapping->m_Dst_from *= 3;
  867.             }
  868.         }
  869.         else if (!width_flag  &&  aln_row.m_Width == 3) {
  870.             // Copy mapping
  871.             mapping.Reset(new CMappingRange(*mappings[map_idx]));
  872.             mapping->m_Src_from *= 3;
  873.             mapping->m_Src_to = (mapping->m_Src_to + 1)*3 - 1;
  874.             mapping->m_Dst_from *= 3;
  875.         }
  877.         if (!mapping->CanMap(start, stop,
  878.             aln_row.m_IsSetStrand, aln_row.m_Strand)) {
  879.             // Mapping does not apply to this segment/row
  880.             continue;
  881.         }
  883.         // Check destination id
  884.         if ( dst_id ) {
  885.             if (mapping->m_Dst_id_Handle != dst_id) {
  886.                 align_flags = eAlign_MultiId;
  887.             }
  888.         }
  889.         dst_id = mapping->m_Dst_id_Handle;
  891.         // At least part of the interval was converted. Calculate
  892.         // trimming coords, trim each row.
  893.         TSeqPos dl = mapping->m_Src_from <= start ?
  894.             0 : mapping->m_Src_from - start;
  895.         TSeqPos dr = mapping->m_Src_to >= stop ?
  896.             0 : stop - mapping->m_Src_to;
  897.         if (dl > 0) {
  898.             // Add segment for the skipped range
  899.             SAlignment_Segment& lseg = x_InsertSeg(seg_it,
  900.                 dl, seg.m_Rows.size());
  901.             for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  902.                 SAlignment_Segment::SAlignment_Row& lrow =
  903.                     lseg.CopyRow(r, seg.m_Rows[r]);
  904.                 if (r == row) {
  905.                     lrow.m_Start = kInvalidSeqPos;
  906.                     lrow.m_Id = dst_id;
  907.                 }
  908.                 else if (lrow.m_Start != kInvalidSeqPos) {
  909.                     int width = seg.m_Rows[r].m_Width ?
  910.                         seg.m_Rows[r].m_Width : 1;
  911.                     if (lrow.SameStrand(aln_row)) {
  912.                         lrow.m_Start += left_shift/width;
  913.                     }
  914.                     else {
  915.                         lrow.m_Start +=
  916.                             (seg.m_Len - lseg.m_Len - left_shift)/width;
  917.                     }
  918.                 }
  919.             }
  920.         }
  921.         start += dl;
  922.         left_shift += dl;
  923.         // At least part of the interval was converted.
  924.         SAlignment_Segment& mseg = x_InsertSeg(seg_it,
  925.             stop - dr - start + 1, seg.m_Rows.size());
  926.         if (!dl  &&  !dr) {
  927.             // copy scores if there's no truncation
  928.             mseg.SetScores(seg.m_Scores);
  929.         }
  930.         for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  931.             SAlignment_Segment::SAlignment_Row& mrow =
  932.                 mseg.CopyRow(r, seg.m_Rows[r]);
  933.             if (r == row) {
  934.                 // translate id and coords
  935.                 CMappingRange::TRange mapped_rg =
  936.                     mapping->Map_Range(start, stop - dr);
  937.                 ENa_strand dst_strand = eNa_strand_unknown;
  938.                 mapping->Map_Strand(
  939.                     aln_row.m_IsSetStrand,
  940.                     aln_row.m_Strand,
  941.                     &dst_strand);
  942.                 mrow.m_Id = mapping->m_Dst_id_Handle;
  943.                 mrow.m_Start = mapped_rg.GetFrom()/dst_width;
  944.                 mrow.m_IsSetStrand |= (dst_strand != eNa_strand_unknown);
  945.                 mrow.m_Strand = dst_strand;
  946.                 mrow.SetMapped();
  947.                 mseg.m_HaveStrands |= mrow.m_IsSetStrand;
  948.             }
  949.             else {
  950.                 if (mrow.m_Start != kInvalidSeqPos) {
  951.                     int width = seg.m_Rows[r].m_Width ?
  952.                         seg.m_Rows[r].m_Width : 1;
  953.                     if (mrow.SameStrand(aln_row)) {
  954.                         mrow.m_Start += left_shift/width;
  955.                     }
  956.                     else {
  957.                         mrow.m_Start +=
  958.                             (seg.m_Len - left_shift - mseg.m_Len)/width;
  959.                     }
  960.                 }
  961.             }
  962.         }
  963.         left_shift += mseg.m_Len;
  964.         start += mseg.m_Len;
  965.         mapped = true;
  966.     }
  967.     if (align_flags == eAlign_MultiId  &&  m_AlignFlags == eAlign_Normal) {
  968.         m_AlignFlags = align_flags;
  969.     }
  970.     if ( !mapped ) {
  971.         // Do not erase the segment, just change the row ID and reset start
  972.         seg.m_Rows[row].m_Start = kInvalidSeqPos;
  973.         seg.m_Rows[row].m_Id = rmap.begin()->second->m_Dst_id_Handle;
  974.         seg.m_Rows[row].SetMapped();
  975.         return seg.m_Rows[row].m_Id;
  976.     }
  977.     if (start <= stop) {
  978.         // Add the remaining unmapped range
  979.         SAlignment_Segment& rseg = x_InsertSeg(seg_it,
  980.             stop - start + 1, seg.m_Rows.size());
  981.         for (size_t r = 0; r < seg.m_Rows.size(); ++r) {
  982.             SAlignment_Segment::SAlignment_Row& rrow =
  983.                 rseg.CopyRow(r, seg.m_Rows[r]);
  984.             if (r == row) {
  985.                 rrow.m_Start = kInvalidSeqPos;
  986.                 rrow.m_Id = dst_id;
  987.             }
  988.             else if (rrow.m_Start != kInvalidSeqPos) {
  989.                 int width = seg.m_Rows[r].m_Width ?
  990.                     seg.m_Rows[r].m_Width : 1;
  991.                 if (rrow.SameStrand(aln_row)) {
  992.                     rrow.m_Start += left_shift/width;
  993.                 }
  994.             }
  995.         }
  996.     }
  997.     m_Segs.erase(old_it);
  998.     return align_flags == eAlign_MultiId ? CSeq_id_Handle() : dst_id;
  999. }
  1002. // Get mapped alignment
  1004. void CSeq_align_Mapper::x_GetDstDendiag(CRef<CSeq_align>& dst) const
  1005. {
  1006.     TDendiag& diags = dst->SetSegs().SetDendiag();
  1007.     ITERATE(TSegments, seg_it, m_Segs) {
  1008.         const SAlignment_Segment& seg = *seg_it;
  1009.         CRef<CDense_diag> diag(new CDense_diag);
  1010.         diag->SetDim(seg.m_Rows.size());
  1011.         diag->SetLen(seg.m_Len);
  1012.         ITERATE(SAlignment_Segment::TRows, row, seg.m_Rows) {
  1013.             CRef<CSeq_id> id(new CSeq_id);
  1014.             id.Reset(&const_cast<CSeq_id&>(*row->m_Id.GetSeqId()));
  1015.             diag->SetIds().push_back(id);
  1016.             diag->SetStarts().push_back(row->GetSegStart());
  1017.             if (seg.m_HaveStrands) { // per-segment strands
  1018.                 diag->SetStrands().push_back(row->m_Strand);
  1019.             }
  1020.         }
  1021.         if ( seg.m_Scores.size() ) {
  1022.             diag->SetScores() = seg.m_Scores;
  1023.         }
  1024.         diags.push_back(diag);
  1025.     }
  1026. }
  1029. void CSeq_align_Mapper::x_GetDstDenseg(CRef<CSeq_align>& dst) const
  1030. {
  1031.     CDense_seg& dseg = dst->SetSegs().SetDenseg();
  1032.     dseg.SetDim(m_Segs.front().m_Rows.size());
  1033.     dseg.SetNumseg(m_Segs.size());
  1034.     if ( m_Segs.front().m_Scores.size() ) {
  1035.         dseg.SetScores() = m_Segs.front().m_Scores;
  1036.     }
  1037.     ITERATE(SAlignment_Segment::TRows, row, m_Segs.front().m_Rows) {
  1038.         CRef<CSeq_id> id(new CSeq_id);
  1039.         id.Reset(&const_cast<CSeq_id&>(*row->m_Id.GetSeqId()));
  1040.         dseg.SetIds().push_back(id);
  1041.         if (m_HaveWidths) {
  1042.             dseg.SetWidths().push_back(row->m_Width);
  1043.         }
  1044.     }
  1045.     ITERATE(TSegments, seg_it, m_Segs) {
  1046.         dseg.SetLens().push_back(seg_it->m_Len);
  1047.         ITERATE(SAlignment_Segment::TRows, row, seg_it->m_Rows) {
  1048.             dseg.SetStarts().push_back(row->GetSegStart());
  1049.             if (m_HaveStrands) { // per-alignment strands
  1050.                 dseg.SetStrands().push_back(row->m_Strand);
  1051.             }
  1052.         }
  1053.     }
  1054. }
  1057. void CSeq_align_Mapper::x_GetDstStd(CRef<CSeq_align>& dst) const
  1058. {
  1059.     TStd& std_segs = dst->SetSegs().SetStd();
  1060.     ITERATE(TSegments, seg_it, m_Segs) {
  1061.         CRef<CStd_seg> std_seg(new CStd_seg);
  1062.         std_seg->SetDim(seg_it->m_Rows.size());
  1063.         if ( seg_it->m_Scores.size() ) {
  1064.             std_seg->SetScores() = seg_it->m_Scores;
  1065.         }
  1066.         ITERATE(SAlignment_Segment::TRows, row, seg_it->m_Rows) {
  1067.             CRef<CSeq_id> id(new CSeq_id);
  1068.             id.Reset(&const_cast<CSeq_id&>(*row->m_Id.GetSeqId()));
  1069.             std_seg->SetIds().push_back(id);
  1070.             CRef<CSeq_loc> loc(new CSeq_loc);
  1071.             if (row->m_Start == kInvalidSeqPos) {
  1072.                 // empty
  1073.                 loc->SetEmpty(*id);
  1074.             }
  1075.             else {
  1076.                 // interval
  1077.                 loc->SetInt().SetId(*id);
  1078.                 TSeqPos len = seg_it->m_Len;
  1079.                 if (row->m_Width) {
  1080.                     len /= row->m_Width;
  1081.                 }
  1082.                 loc->SetInt().SetFrom(row->m_Start);
  1083.                 // len may be 0 after dividing by width
  1084.                 loc->SetInt().SetTo(row->m_Start + (len ? len - 1 : 0));
  1085.                 if (row->m_IsSetStrand) {
  1086.                     loc->SetInt().SetStrand(row->m_Strand);
  1087.                 }
  1088.             }
  1089.             std_seg->SetLoc().push_back(loc);
  1090.         }
  1091.         std_segs.push_back(std_seg);
  1092.     }
  1093. }
  1096. void CSeq_align_Mapper::x_GetDstPacked(CRef<CSeq_align>& dst) const
  1097. {
  1098.     CPacked_seg& pseg = dst->SetSegs().SetPacked();
  1099.     pseg.SetDim(m_Segs.front().m_Rows.size());
  1100.     pseg.SetNumseg(m_Segs.size());
  1101.     if ( m_Segs.front().m_Scores.size() ) {
  1102.         pseg.SetScores() = m_Segs.front().m_Scores;
  1103.     }
  1104.     ITERATE(SAlignment_Segment::TRows, row, m_Segs.front().m_Rows) {
  1105.         CRef<CSeq_id> id(new CSeq_id);
  1106.         id.Reset(&const_cast<CSeq_id&>(*row->m_Id.GetSeqId()));
  1107.         pseg.SetIds().push_back(id);
  1108.     }
  1109.     ITERATE(TSegments, seg_it, m_Segs) {
  1110.         pseg.SetLens().push_back(seg_it->m_Len);
  1111.         ITERATE(SAlignment_Segment::TRows, row, seg_it->m_Rows) {
  1112.             pseg.SetStarts().push_back(row->GetSegStart());
  1113.             pseg.SetPresent().push_back(row->m_Start != kInvalidSeqPos);
  1114.             if (m_HaveStrands) {
  1115.                 pseg.SetStrands().push_back(row->m_Strand);
  1116.             }
  1117.         }
  1118.     }
  1119. }
  1122. void CSeq_align_Mapper::x_GetDstDisc(CRef<CSeq_align>& dst) const
  1123. {
  1124.     CSeq_align_set::Tdata& data = dst->SetSegs().SetDisc().Set();
  1125.     ITERATE(TSubAligns, it, m_SubAligns) {
  1126.         data.push_back(it->GetDstAlign());
  1127.     }
  1128. }
  1131. CRef<CSeq_align> CSeq_align_Mapper::GetDstAlign(void) const
  1132. {
  1133.     if (m_DstAlign) {
  1134.         return m_DstAlign;
  1135.     }
  1137.     CRef<CSeq_align> dst(new CSeq_align);
  1138.     dst->SetType(m_OrigAlign->GetType());
  1139.     if (m_OrigAlign->IsSetDim()) {
  1140.         dst->SetDim(m_OrigAlign->GetDim());
  1141.     }
  1142.     if (m_OrigAlign->IsSetScore()) {
  1143.         dst->SetScore() = m_OrigAlign->GetScore();
  1144.     }
  1145.     if (m_OrigAlign->IsSetBounds()) {
  1146.         dst->SetBounds() = m_OrigAlign->GetBounds();
  1147.     }
  1148.     switch ( m_OrigAlign->GetSegs().Which() ) {
  1149.     case CSeq_align::C_Segs::e_Dendiag:
  1150.         {
  1151.             x_GetDstDendiag(dst);
  1152.             break;
  1153.         }
  1154.     case CSeq_align::C_Segs::e_Denseg:
  1155.         {
  1156.             if (m_AlignFlags == eAlign_Normal) {
  1157.                 x_GetDstDenseg(dst);
  1158.             }
  1159.             else {
  1160.                 x_GetDstDendiag(dst);
  1161.             }
  1162.             break;
  1163.         }
  1164.     case CSeq_align::C_Segs::e_Std:
  1165.         {
  1166.             x_GetDstStd(dst);
  1167.             break;
  1168.         }
  1169.     case CSeq_align::C_Segs::e_Packed:
  1170.         {
  1171.             if (m_AlignFlags == eAlign_Normal) {
  1172.                 x_GetDstPacked(dst);
  1173.             }
  1174.             else {
  1175.                 x_GetDstDendiag(dst);
  1176.             }
  1177.             break;
  1178.         }
  1179.     case CSeq_align::C_Segs::e_Disc:
  1180.         {
  1181.             x_GetDstDisc(dst);
  1182.             break;
  1183.         }
  1184.     default:
  1185.         {
  1186.             dst->Assign(*m_OrigAlign);
  1187.             break;
  1188.         }
  1189.     }
  1190.     return m_DstAlign = dst;
  1191. }
  1194. END_SCOPE(objects)
  1195. END_NCBI_SCOPE
  1197. /*
  1198. * ---------------------------------------------------------------------------
  1199. * $Log: seq_align_mapper.cpp,v $
  1200. * Revision 1000.1  2004/06/01 19:23:43  gouriano
  1201. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.10
  1202. *
  1203. * Revision 1.10  2004/05/27 17:52:11  grichenk
  1204. * Fixed warnings
  1205. *
  1206. * Revision 1.9  2004/05/26 14:57:47  ucko
  1207. * Change some types from unsigned int to size_t for consistency with STL
  1208. * containers' size() methods.
  1209. *
  1210. * Revision 1.8  2004/05/26 14:29:20  grichenk
  1211. * Redesigned CSeq_align_Mapper: preserve non-mapping intervals,
  1212. * fixed strands handling, improved performance.
  1213. *
  1214. * Revision 1.7  2004/05/21 21:42:12  gorelenk
  1215. * Added PCH ncbi_pch.hpp
  1216. *
  1217. * Revision 1.6  2004/05/13 19:10:57  grichenk
  1218. * Preserve scores in mapped alignments
  1219. *
  1220. * Revision 1.5  2004/05/10 20:22:24  grichenk
  1221. * Fixed more warnings (removed inlines)
  1222. *
  1223. * Revision 1.4  2004/04/12 14:35:59  grichenk
  1224. * Fixed mapping of alignments between nucleotides and proteins
  1225. *
  1226. * Revision 1.3  2004/04/07 18:36:12  grichenk
  1227. * Fixed std-seg mapping
  1228. *
  1229. * Revision 1.2  2004/04/01 20:11:17  rsmith
  1230. * add missing break to switch on seq-id types.
  1231. *
  1232. * Revision 1.1  2004/03/30 15:40:35  grichenk
  1233. * Initial revision
  1234. *
  1235. *
  1236. * ===========================================================================
  1237. */