seq_loc_mapper.cpp

来自「ncbi源码」· C++ 代码 · 共 1,895 行 · 第 1/5 页

/*
 * ===========================================================================
 * PRODUCTION $Log: seq_loc_mapper.cpp,v $
 * PRODUCTION Revision 1000.1  2004/06/01 19:24:12  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.20
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: seq_loc_mapper.cpp,v 1000.1 2004/06/01 19:24:12 gouriano Exp $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's official duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================
*
* Author: Aleksey Grichenko
*
* File Description:
*   Seq-loc mapper
*
*/

#include <ncbi_pch.hpp>
#include <objmgr/seq_loc_mapper.hpp>
#include <objmgr/scope.hpp>
#include <objmgr/objmgr_exception.hpp>
#include <objmgr/seq_map.hpp>
#include <objmgr/seq_map_ci.hpp>
#include <objmgr/impl/synonyms.hpp>
#include <objmgr/impl/seq_align_mapper.hpp>
#include <objmgr/impl/seq_loc_cvt.hpp>
#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqfeat/Seq_feat.hpp>
#include <objects/seqfeat/Cdregion.hpp>
#include <objects/seqloc/Seq_loc_equiv.hpp>
#include <objects/seqloc/Seq_bond.hpp>
#include <objects/seqalign/seqalign__.hpp>
#include <algorithm>


BEGIN_NCBI_SCOPE
BEGIN_SCOPE(objects)


CMappingRange::CMappingRange(CSeq_id_Handle     src_id,
                             TSeqPos            src_from,
                             TSeqPos            src_length,
                             ENa_strand         src_strand,
                             CSeq_id_Handle     dst_id,
                             TSeqPos            dst_from,
                             ENa_strand         dst_strand)
    : m_Src_id_Handle(src_id),
      m_Src_from(src_from),
      m_Src_to(src_from + src_length - 1),
      m_Src_strand(src_strand),
      m_Dst_id_Handle(dst_id),
      m_Dst_from(dst_from),
      m_Dst_strand(dst_strand),
      m_Reverse(!SameOrientation(src_strand, dst_strand))
{
    return;
}


bool CMappingRange::CanMap(TSeqPos from,
                           TSeqPos to,
                           bool is_set_strand,
                           ENa_strand strand) const
{
    if (from > m_Src_to || to < m_Src_from) {
        return false;
    }
    return SameOrientation(is_set_strand ? strand : eNa_strand_unknown,
        m_Src_strand)  ||  (m_Src_strand == eNa_strand_unknown);
}


TSeqPos CMappingRange::Map_Pos(TSeqPos pos) const
{
    _ASSERT(pos >= m_Src_from  &&  pos <= m_Src_to);
    if (!m_Reverse) {
        return m_Dst_from + pos - m_Src_from;
    }
    else {
        return m_Dst_from + m_Src_to - pos;
    }
}


CMappingRange::TRange CMappingRange::Map_Range(TSeqPos from, TSeqPos to) const
{
    if (!m_Reverse) {
        return TRange(Map_Pos(max(from, m_Src_from)),
            Map_Pos(min(to, m_Src_to)));
    }
    else {
        return TRange(Map_Pos(min(to, m_Src_to)),
            Map_Pos(max(from, m_Src_from)));
    }
}


bool CMappingRange::Map_Strand(bool is_set_strand,
                               ENa_strand src,
                               ENa_strand* dst) const
{
    _ASSERT(dst);
    if (is_set_strand) {
        *dst = m_Reverse ? Reverse(src) : src;
        return true;
    }
    if (m_Dst_strand != eNa_strand_unknown) {
        // Destination strand may be set for nucleotides
        // even if the source one is not set.
        *dst = m_Dst_strand;
        return true;
    }
    return false;
}


const CMappingRange::TFuzz kEmptyFuzz(0);

CInt_fuzz::ELim CMappingRange::x_ReverseFuzzLim(CInt_fuzz::ELim lim) const
{
    switch ( lim ) {
    case CInt_fuzz::eLim_gt:
        return CInt_fuzz::eLim_lt;
    case CInt_fuzz::eLim_lt:
        return CInt_fuzz::eLim_gt;
    case CInt_fuzz::eLim_tr:
        return CInt_fuzz::eLim_tl;
    case CInt_fuzz::eLim_tl:
        return CInt_fuzz::eLim_tr;
    default:
        return lim;
    }
}


CMappingRange::TRangeFuzz CMappingRange::Map_Fuzz(TRangeFuzz& fuzz) const
{
    // Maps some fuzz types to reverse strand
    if ( !m_Reverse ) {
        return fuzz;
    }
    TRangeFuzz res(fuzz.second, fuzz.first);
    if ( res.first ) {
        switch ( res.first->Which() ) {
        case CInt_fuzz::e_Lim:
            {
                res.first->SetLim(x_ReverseFuzzLim(res.first->GetLim()));
                break;
            }
        default:
            // Other types are not converted
            break;
        }
    }
    if ( fuzz.second ) {
        switch ( res.second->Which() ) {
        case CInt_fuzz::e_Lim:
            {
                res.second->SetLim(x_ReverseFuzzLim(res.second->GetLim()));
                break;
            }
        default:
            // Other types are not converted
            break;
        }
    }
    return res;
}


CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_feat&  map_feat,
                                 EFeatMapDirection dir,
                                 CScope*           scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    _ASSERT(map_feat.IsSetProduct());
    int frame = 0;
    if (map_feat.GetData().IsCdregion()) {
        frame = map_feat.GetData().GetCdregion().GetFrame();
    }
    if (dir == eLocationToProduct) {
        x_Initialize(map_feat.GetLocation(), map_feat.GetProduct(), frame);
    }
    else {
        x_Initialize(map_feat.GetProduct(), map_feat.GetLocation(), frame);
    }
}


CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_loc& source,
                                 const CSeq_loc& target,
                                 CScope* scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    x_Initialize(source, target);
}


CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_align& map_align,
                                 const CSeq_id&    to_id,
                                 CScope*           scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    x_Initialize(map_align, to_id);
}


CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeq_align& map_align,
                                 size_t            to_row,
                                 CScope*           scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    x_Initialize(map_align, to_row);
}


CSeq_loc_Mapper::CSeq_loc_Mapper(CBioseq_Handle target_seq)
    : m_Scope(&target_seq.GetScope()),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    CConstRef<CSeq_id> dst_id = target_seq.GetSeqId();
    x_Initialize(target_seq.GetSeqMap(),
        dst_id.GetPointerOrNull());
    // Ignore seq-map destination ranges, map whole sequence to itself,
    // use unknown strand only.
    m_DstRanges.resize(1);
    m_DstRanges[0].clear();
    if (m_Scope) {
        CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*dst_id);
        ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
            m_DstRanges[0][target_seq.GetSeq_id_Handle()]
                .push_back(TRange::GetWhole());
        }
    }
    else {
        m_DstRanges[0][target_seq.GetSeq_id_Handle()]
            .push_back(TRange::GetWhole());
    }
}


CSeq_loc_Mapper::CSeq_loc_Mapper(const CSeqMap& seq_map,
                                 const CSeq_id* dst_id,
                                 CScope*        scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    x_Initialize(seq_map, dst_id);
}


CSeq_loc_Mapper::CSeq_loc_Mapper(size_t          depth,
                                 CBioseq_Handle& source_seq)
    : m_Scope(&source_seq.GetScope()),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    if (depth > 0) {
        depth--;
        x_Initialize(source_seq.GetSeqMap(),
            depth, source_seq.GetSeqId().GetPointer());
    }
    else /* if (depth == 0) */ {
        // Synonyms conversion
        CConstRef<CSeq_id> dst_id = source_seq.GetSeqId();
        m_DstRanges.resize(1);
        if (m_Scope) {
            CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*dst_id);
            ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
                m_DstRanges[0][source_seq.GetSeq_id_Handle()]
                    .push_back(TRange::GetWhole());
            }
        }
        else {
            m_DstRanges[0][source_seq.GetSeq_id_Handle()]
                .push_back(TRange::GetWhole());
        }
    }
}


CSeq_loc_Mapper::CSeq_loc_Mapper(size_t         depth,
                                 const CSeqMap& source_seqmap,
                                 const CSeq_id* src_id,
                                 CScope*        scope)
    : m_Scope(scope),
      m_MergeFlag(eMergeNone),
      m_GapFlag(eGapPreserve),
      m_KeepNonmapping(false),
      m_UseWidth(false),
      m_Dst_width(0)
{
    if (depth > 0) {
        depth--;
        x_Initialize(source_seqmap, depth, src_id);
    }
    else /* if (depth == 0) */ {
        // Synonyms conversion
        m_DstRanges.resize(1);
        if (bool(m_Scope)  &&  src_id) {
            CSeq_id_Handle src_idh = CSeq_id_Handle::GetHandle(*src_id);
            CConstRef<CSynonymsSet> dst_syns = m_Scope->GetSynonyms(*src_id);
            ITERATE(CSynonymsSet, dst_syn_it, *dst_syns) {
                m_DstRanges[0][src_idh]
                    .push_back(TRange::GetWhole());
            }
        }
    }
}


CSeq_loc_Mapper::~CSeq_loc_Mapper(void)
{
    return;
}


void CSeq_loc_Mapper::PreserveDestinationLocs(void)
{
    for (size_t str_idx = 0; str_idx < m_DstRanges.size(); str_idx++) {
        NON_CONST_ITERATE(TDstIdMap, id_it, m_DstRanges[str_idx]) {
            id_it->second.sort();
            TSeqPos dst_start = kInvalidSeqPos;
            TSeqPos dst_stop = kInvalidSeqPos;