annot_collector.cpp

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

    // gather x annotation type
    const CAnnotObject_Info* x_info;
    CSeq_annot::C_Data::E_Choice x_annot_type;
    if ( !x.IsSNPFeat() ) {
        const CSeq_annot_Info& x_annot = x.GetSeq_annot_Info();
        x_info = &x_annot.GetAnnotObject_Info(x.GetAnnotObjectIndex());
        x_annot_type = x_info->GetAnnotType();
    }
    else {
        x_info = 0;
        x_annot_type = CSeq_annot::C_Data::e_Ftable;
    }

    // gather y annotation type
    const CAnnotObject_Info* y_info;
    CSeq_annot::C_Data::E_Choice y_annot_type;
    if ( !y.IsSNPFeat() ) {
        const CSeq_annot_Info& y_annot = y.GetSeq_annot_Info();
        y_info = &y_annot.GetAnnotObject_Info(y.GetAnnotObjectIndex());
        y_annot_type = y_info->GetAnnotType();
    }
    else {
        y_info = 0;
        y_annot_type = CSeq_annot::C_Data::e_Ftable;
    }

    // compare by annotation type (feature, align, graph)
    if ( x_annot_type != y_annot_type ) {
        return x_annot_type < y_annot_type;
    }

    if ( x_annot_type == CSeq_annot::C_Data::e_Ftable ) {
        // compare features by type
        if ( !x_info != !y_info ) {
            CSeqFeatData::E_Choice x_feat_type = CSeqFeatData::e_Imp;
            CSeqFeatData::ESubtype x_feat_subtype =
                CSeqFeatData::eSubtype_variation;
            if ( x_info ) {
                x_feat_type = x_info->GetFeatType();
                x_feat_subtype = x_info->GetFeatSubtype();
            }

            CSeqFeatData::E_Choice y_feat_type = CSeqFeatData::e_Imp;
            CSeqFeatData::ESubtype y_feat_subtype =
                CSeqFeatData::eSubtype_variation;
            if ( y_info ) {
                y_feat_type = y_info->GetFeatType();
                y_feat_subtype = y_info->GetFeatSubtype();
            }

            // one is simple SNP feature, another is some complex feature
            if ( x_feat_type != y_feat_type ) {
                int x_order = CSeq_feat::GetTypeSortingOrder(x_feat_type);
                int y_order = CSeq_feat::GetTypeSortingOrder(y_feat_type);
                if ( x_order != y_order ) {
                    return x_order < y_order;
                }
            }

            // non-variations first
            if ( x_feat_subtype != y_feat_subtype ) {
                return x_feat_subtype < y_feat_subtype;
            }

            // both are variations but one is simple and another is not.
            // required order: simple first.
            // if !x_info == true -> x - simple, y - complex -> return true
            // if !x_info == false -> x - complex, y - simple -> return false
            return !x_info;
        }

        if ( x_info ) {
            // both are complex features
            try {
                const CSeq_feat* x_feat = x_info->GetFeatFast();
                const CSeq_feat* y_feat = y_info->GetFeatFast();
                _ASSERT(x_feat && y_feat);
                int diff = x_feat->CompareNonLocation(*y_feat,
                                                      GetLocation(x, *x_feat),
                                                      GetLocation(y, *y_feat));
                if ( diff != 0 ) {
                    return diff < 0;
                }
            }
            catch ( exception& /*ignored*/ ) {
                // do not fail sort when compare function throws an exception
            }
        }
    }
    return x < y;
}


struct CAnnotObject_Less
{
    // Compare CRef-s: both must be features
    bool operator()(const CAnnotObject_Ref& x,
                    const CAnnotObject_Ref& y) const
        {
            {
                // smallest left extreme first
                TSeqPos x_from = x.GetFrom();
                TSeqPos y_from = y.GetFrom();
                if ( x_from != y_from ) {
                    return x_from < y_from;
                }
            }
            {
                // longest feature first
                TSeqPos x_to = x.GetToOpen();
                TSeqPos y_to = y.GetToOpen();
                if ( x_to != y_to ) {
                    return x_to > y_to;
                }
            }
            return type_less(x, y);
        }
    CAnnotObjectType_Less type_less;
};


struct CAnnotObject_LessReverse
{
    // Compare CRef-s: both must be features
    bool operator()(const CAnnotObject_Ref& x,
                    const CAnnotObject_Ref& y) const
        {
            {
                // largest right extreme first
                TSeqPos x_to = x.GetToOpen();
                TSeqPos y_to = y.GetToOpen();
                if ( x_to != y_to ) {
                    return x_to > y_to;
                }
            }
            {
                // longest feature first
                TSeqPos x_from = x.GetFrom();
                TSeqPos y_from = y.GetFrom();
                if ( x_from != y_from ) {
                    return x_from < y_from;
                }
            }
            return type_less(x, y);
        }
    CAnnotObjectType_Less type_less;
};


/////////////////////////////////////////////////////////////////////////////
// CAnnot_Collector, CAnnotMappingCollector
/////////////////////////////////////////////////////////////////////////////


class CAnnotMappingCollector
{
public:
    typedef map<CAnnotObject_Ref,
                CRef<CSeq_loc_Conversion_Set> > TAnnotMappingSet;
    // Set of annotations for complex remapping
    TAnnotMappingSet              m_AnnotMappingSet;
    // info of limit object
    CConstRef<CObject>            m_LimitObjectInfo;
};


CAnnot_Collector::CAnnot_Collector(const SAnnotSelector& selector,
                                         CScope&               scope)
    : m_Selector(selector),
      m_Scope(scope),
      m_MappingCollector(new CAnnotMappingCollector)
{
    return;
}

CAnnot_Collector::~CAnnot_Collector(void)
{
    x_Clear();
}


inline
size_t CAnnot_Collector::x_GetAnnotCount(void) const
{
    return m_AnnotSet.size() +
        (m_MappingCollector.get() ?
        m_MappingCollector->m_AnnotMappingSet.size() : 0);
}


void CAnnot_Collector::x_Clear(void)
{
    m_AnnotSet.clear();
    if ( m_MappingCollector.get() ) {
        m_MappingCollector.reset();
    }
    m_TSE_LockSet.clear();
    m_Scope = CHeapScope();
}


void CAnnot_Collector::x_Initialize(const CBioseq_Handle& bioseq,
                                       TSeqPos start, TSeqPos stop)
{
    try {
        if ( start == 0 && stop == 0 ) {
            stop = bioseq.GetBioseqLength()-1;
        }
        CHandleRangeMap master_loc;
        master_loc.AddRange(bioseq.GetSeq_id_Handle(),
                            CHandleRange::TRange(start, stop),
                            eNa_strand_unknown);
        x_Initialize(master_loc);
    } catch (...) {
        // clear all members - GCC 3.0.4 does not do it
        x_Clear();
        throw;
    }
}


void CAnnot_Collector::x_Initialize(const CHandleRangeMap& master_loc)
{
    try {
        if ( !m_Selector.m_LimitObject ) {
            m_Selector.m_LimitObjectType = SAnnotSelector::eLimit_None;
        }
        if ( m_Selector.m_LimitObjectType != SAnnotSelector::eLimit_None ) {
            x_GetTSE_Info();
        }

        bool found = x_Search(master_loc, 0);
        bool deeper = !(found && m_Selector.m_AdaptiveDepth) &&
            m_Selector.m_ResolveMethod != SAnnotSelector::eResolve_None  &&
            m_Selector.m_ResolveDepth > 0;
        if ( deeper ) {
            ITERATE ( CHandleRangeMap::TLocMap, idit, master_loc.GetMap() ) {
                //### Check for eLoadedOnly
                CBioseq_Handle bh = m_Scope->GetBioseqHandle(idit->first,
                    CScope::eGetBioseq_All);
                if ( !bh ) {
                    if (m_Selector.m_IdResolving ==
                        SAnnotSelector::eFailUnresolved) {
                        // resolve by Seq-id only
                        NCBI_THROW(CAnnotException, eFindFailed,
                                   "Cannot resolve master id");
                    }
                    continue; // Do not crash - just skip unresolvable IDs
                }

                CRef<CSeq_loc> master_loc_empty(new CSeq_loc);
                master_loc_empty->SetEmpty(
                    const_cast<CSeq_id&>(*idit->first.GetSeqId()));
                
                CHandleRange::TRange idrange =
                    idit->second.GetOverlappingRange();
                const CSeqMap& seqMap = bh.GetSeqMap();
                CSeqMap_CI smit(seqMap.FindResolved(idrange.GetFrom(),
                                                    m_Scope,
                                                    m_Selector.m_ResolveDepth-1,
                                                    CSeqMap::fFindRef));
                while ( smit && smit.GetPosition() < idrange.GetToOpen() ) {
                    _ASSERT(smit.GetType() == CSeqMap::eSeqRef);
                    if ( m_Selector.m_ResolveMethod ==
                        SAnnotSelector::eResolve_TSE &&
                        !m_Scope->GetBioseqHandleFromTSE(smit.GetRefSeqid(),
                                                         bh) ) {
                        smit.Next(false);
                        continue;
                    }
                    found = x_SearchMapped(smit,
                                           *master_loc_empty,
                                           idit->first,
                                           idit->second);
                    deeper = !(found && m_Selector.m_AdaptiveDepth);
                    smit.Next(deeper);
                }
            }
        }
        NON_CONST_ITERATE(CAnnotMappingCollector::TAnnotMappingSet, amit,
            m_MappingCollector->m_AnnotMappingSet) {
            CAnnotObject_Ref annot_ref = amit->first;
            amit->second->Convert(annot_ref,
                m_Selector.m_FeatProduct ? CSeq_loc_Conversion::eProduct :
                CSeq_loc_Conversion::eLocation);
            m_AnnotSet.push_back(annot_ref);
        }
        m_MappingCollector->m_AnnotMappingSet.clear();
        x_Sort();
        m_MappingCollector.reset();
    }
    catch (...) {
        // clear all members - GCC 3.0.4 does not do it
        x_Clear();
        throw;
    }
}


void CAnnot_Collector::x_Initialize(void)
{
    try {
        // Limit must be set, resolving is obsolete
        _ASSERT(m_Selector.m_LimitObjectType != SAnnotSelector::eLimit_None);
        _ASSERT(m_Selector.m_LimitObject);
        _ASSERT(m_Selector.m_ResolveMethod == SAnnotSelector::eResolve_None);
        x_GetTSE_Info();

        x_SearchAll();
        x_Sort();
        m_MappingCollector.reset();
    }
    catch (...) {
        // clear all members - GCC 3.0.4 does not do it
        x_Clear();
        throw;
    }
}


void CAnnot_Collector::x_Sort(void)
{
    _ASSERT(m_MappingCollector->m_AnnotMappingSet.empty());
    switch ( m_Selector.m_SortOrder ) {
    case SAnnotSelector::eSortOrder_Normal:
        sort(m_AnnotSet.begin(), m_AnnotSet.end(), CAnnotObject_Less());
        break;
    case SAnnotSelector::eSortOrder_Reverse:
        sort(m_AnnotSet.begin(), m_AnnotSet.end(), CAnnotObject_LessReverse());
        break;
    default:
        // do nothing
        break;
    }
}


inline
bool
CAnnot_Collector::x_MatchLimitObject(const CAnnotObject_Info& object_info) const
{
    if ( m_Selector.m_LimitObjectType == SAnnotSelector::eLimit_Seq_entry_Info ) {
        const CSeq_entry_Info* info = &object_info.GetSeq_entry_Info();
        _ASSERT(m_MappingCollector->m_LimitObjectInfo);
        _ASSERT(info);
        for ( ;; ) {
            if ( info == m_MappingCollector->m_LimitObjectInfo.GetPointer() ) {
                return true;
            }
            if ( !info->HasParent_Info() ) {
                return false;
            }
            info = &info->GetParentSeq_entry_Info();
        }
    }
    else if ( m_Selector.m_LimitObjectType ==
        SAnnotSelector::eLimit_Seq_annot_Info ) {
        const CSeq_annot_Info* info = &object_info.GetSeq_annot_Info();