seq_map.cpp

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

/*
 * ===========================================================================
 * PRODUCTION $Log: seq_map.cpp,v $
 * PRODUCTION Revision 1000.4  2004/06/01 19:24:15  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.53
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: seq_map.cpp,v 1000.4 2004/06/01 19:24:15 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.
*
* ===========================================================================
*
* Authors: Aleksey Grichenko, Michael Kimelman, Eugene Vasilchenko,
*          Andrei Gourianov
*
* File Description:
*   Sequence map for the Object Manager. Describes sequence as a set of
*   segments of different types (data, reference, gap or end).
*
*/

#include <ncbi_pch.hpp>
#include <objmgr/seq_map.hpp>
#include <objmgr/seq_map_ci.hpp>
#include <objmgr/seq_map_ext.hpp>
#include <objmgr/seq_id_mapper.hpp>
#include <objmgr/scope.hpp>
#include <objmgr/bioseq_handle.hpp>

#include <objects/seq/Bioseq.hpp>
#include <objects/seq/Seq_data.hpp>
#include <objects/seq/Seq_inst.hpp>
#include <objects/seq/Delta_ext.hpp>
#include <objects/seq/Delta_seq.hpp>
#include <objects/seq/Seq_literal.hpp>
#include <objects/seq/Seq_ext.hpp>
#include <objects/seq/Seg_ext.hpp>
#include <objects/seq/Ref_ext.hpp>

#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqloc/Seq_point.hpp>
#include <objects/seqloc/Seq_loc_mix.hpp>
#include <objects/seqloc/Seq_loc_equiv.hpp>
#include <objects/seqloc/Seq_interval.hpp>
#include <objects/seqloc/Packed_seqint.hpp>
#include <objects/seqloc/Packed_seqpnt.hpp>

#include <algorithm>

BEGIN_NCBI_SCOPE
BEGIN_SCOPE(objects)

////////////////////////////////////////////////////////////////////
//  CSeqMap::CSegment

inline
CSeqMap::CSegment::CSegment(ESegmentType seg_type,
                            TSeqPos length)
    : m_Position(kInvalidSeqPos),
      m_Length(length),
      m_SegType(seg_type),
      m_RefMinusStrand(false),
      m_RefPosition(0)
{
}

////////////////////////////////////////////////////////////////////
//  CSeqMap


CSeqMap::CSeqMap(void)
    : m_Resolved(0),
      m_Mol(CSeq_inst::eMol_not_set),
      m_SeqLength(kInvalidSeqPos)
{
}


CSeqMap::CSeqMap(CSeqMap* /*parent*/, size_t /*index*/)
    : m_Resolved(0),
      m_Mol(CSeq_inst::eMol_not_set),
      m_SeqLength(kInvalidSeqPos)
{
}


CSeqMap::CSeqMap(const CSeq_loc& ref)
    : m_Resolved(0),
      m_Mol(CSeq_inst::eMol_not_set),
      m_SeqLength(kInvalidSeqPos)
{
    x_AddEnd();
    x_Add(ref);
    x_AddEnd();
}


CSeqMap::CSeqMap(const CSeq_data& data, TSeqPos length)
    : m_Resolved(0),
      m_Mol(CSeq_inst::eMol_not_set),
      m_SeqLength(kInvalidSeqPos)
{
    x_AddEnd();
    x_Add(data, length);
    x_AddEnd();
}


CSeqMap::CSeqMap(TSeqPos length)
    : m_Resolved(0),
      m_Mol(CSeq_inst::eMol_not_set),
      m_SeqLength(length)
{
    x_AddEnd();
    x_AddGap(length);
    x_AddEnd();
}


CSeqMap::~CSeqMap(void)
{
    m_Resolved = 0;
    m_Segments.clear();
}


void CSeqMap::x_GetSegmentException(size_t /*index*/) const
{
    NCBI_THROW(CSeqMapException, eInvalidIndex,
               "Invalid segment index");
}


CSeqMap::CSegment& CSeqMap::x_SetSegment(size_t index)
{
    if ( index >= x_GetSegmentsCount() ) {
        NCBI_THROW(CSeqMapException, eInvalidIndex,
                   "Invalid segment index");
    }
    return m_Segments[index];
}


CBioseq_Handle CSeqMap::x_GetBioseqHandle(const CSegment& seg, CScope* scope) const
{
    if ( !scope ) {
        NCBI_THROW(CSeqMapException, eNullPointer,
                   "Null scope pointer");
    }
    CBioseq_Handle bh = scope->GetBioseqHandle(x_GetRefSeqid(seg));
    if ( !bh ) {
        NCBI_THROW(CSeqMapException, eFail,
                   "Cannot resolve reference");
    }
    return bh;
}


TSeqPos CSeqMap::x_ResolveSegmentLength(size_t index, CScope* scope) const
{
    const CSegment& seg = x_GetSegment(index);
    TSeqPos length = seg.m_Length;
    if ( length == kInvalidSeqPos ) {
        if ( seg.m_SegType == eSeqSubMap ) {
            length = x_GetSubSeqMap(seg, scope)->GetLength(scope);
        }
        else if ( seg.m_SegType == eSeqRef ) {
            length = x_GetBioseqHandle(seg, scope).GetBioseqLength();
        }
        _ASSERT(length != kInvalidSeqPos);
        seg.m_Length = length;
    }
    return length;
}


TSeqPos CSeqMap::x_ResolveSegmentPosition(size_t index, CScope* scope) const
{
    if ( index > x_GetSegmentsCount() ) {
        NCBI_THROW(CSeqMapException, eInvalidIndex,
                   "Invalid segment index");
    }
    size_t resolved = m_Resolved;
    if ( index <= resolved )
        return x_GetSegment(index).m_Position;
    TSeqPos resolved_pos = x_GetSegment(resolved).m_Position;
    do {
        resolved_pos += x_GetSegmentLength(resolved, scope);
        m_Segments[++resolved].m_Position = resolved_pos;
    } while ( resolved < index );
    {{
        CFastMutexGuard guard(m_SeqMap_Mtx);
        if ( m_Resolved < resolved )
            m_Resolved = resolved;
    }}
    return resolved_pos;
}


size_t CSeqMap::x_FindSegment(TSeqPos pos, CScope* scope) const
{
    size_t resolved = m_Resolved;
    TSeqPos resolved_pos = x_GetSegment(resolved).m_Position;
    if ( resolved_pos <= pos ) {
        do {
            if ( resolved >= x_GetSegmentsCount() ) {
                // end of segments
                m_Resolved = resolved;
                return size_t(-1);
            }
            resolved_pos += x_GetSegmentLength(resolved, scope);
            m_Segments[++resolved].m_Position = resolved_pos;
        } while ( resolved_pos <= pos );
        {{
            CFastMutexGuard guard(m_SeqMap_Mtx);
            if ( m_Resolved < resolved )
                m_Resolved = resolved;
        }}
        return resolved - 1;
    }
    else {
        TSegments::const_iterator end = m_Segments.begin()+resolved;
        TSegments::const_iterator it = 
            upper_bound(m_Segments.begin(), end,
                        pos, SPosLessSegment());
        if ( it == end ) {
            return size_t(-1);
        }
        return it - m_Segments.begin();
    }
}


void CSeqMap::x_LoadObject(const CSegment& seg) const
{
    _ASSERT(seg.m_Position != kInvalidSeqPos);
    if ( !seg.m_RefObject ) {
        NCBI_THROW(CSeqMapException, eFail, "Cannot load part of seqmap");
    }
}


CConstRef<CSeqMap> CSeqMap::x_GetSubSeqMap(const CSegment& seg, CScope* scope,
                                           bool resolveExternal) const
{
    CConstRef<CSeqMap> ret;
    if ( seg.m_SegType == eSeqSubMap ) {
        if ( !seg.m_RefObject ) {
            x_LoadObject(seg);
        }
        if ( !seg.m_RefObject ) {
        NCBI_THROW(CSeqMapException, eNullPointer, "Null CSeqMap pointer");
        }
        ret.Reset(static_cast<const CSeqMap*>(seg.m_RefObject.GetPointer()));
    }
    else if ( resolveExternal && seg.m_SegType == eSeqRef ) {
        ret.Reset(&x_GetBioseqHandle(seg, scope).GetSeqMap());
    }
    return ret;
}


void CSeqMap::x_SetSubSeqMap(size_t /*index*/, CSeqMap_Delta_seqs* /*subMap*/)
{
    // not valid in generic seq map -> incompatible objects
    NCBI_THROW(CSeqMapException, eDataError, "Invalid parent map");
}


const CSeq_data& CSeqMap::x_GetSeq_data(const CSegment& seg) const
{
    if ( seg.m_SegType == eSeqData ) {
        if ( !seg.m_RefObject ) {
            x_LoadObject(seg);
        }
        if ( !seg.m_RefObject ) {
            NCBI_THROW(CSeqMapException, eNullPointer,
                       "Null CSeq_data pointer");
        }
        return *static_cast<const CSeq_data*>(seg.m_RefObject.GetPointer());
    }
    NCBI_THROW(CSeqMapException, eSegmentTypeError,
               "Invalid segment type");
}


void CSeqMap::x_SetSeq_data(size_t index, CSeq_data& data)
{
    // check segment type
    CSegment& seg = x_SetSegment(index);
    if ( seg.m_SegType != eSeqData ) {
        NCBI_THROW(CSeqMapException, eSegmentTypeError,
                   "Invalid segment type");
    }
    // lock for object modification
    CFastMutexGuard guard(m_SeqMap_Mtx);
    // check for object
    if ( seg.m_RefObject ) {
        NCBI_THROW(CSeqMapException, eDataError, "CSeq_data already set");
    }
    // set object
    seg.m_RefObject.Reset(&data);
}


const CSeq_id& CSeqMap::x_GetRefSeqid(const CSegment& seg) const
{
    if ( seg.m_SegType == eSeqRef ) {
        if ( !seg.m_RefObject ) {
            NCBI_THROW(CSeqMapException, eNullPointer, "Null CSeq_id pointer");
        }
        return static_cast<const CSeq_id&>(*seg.m_RefObject);
    }
    NCBI_THROW(CSeqMapException, eSegmentTypeError,
               "Invalid segment type");
}


TSeqPos CSeqMap::x_GetRefPosition(const CSegment& seg) const
{
    return seg.m_RefPosition;