density_map.hpp

来自「ncbi源码」· HPP 代码 · 共 662 行 · 第 1/2 页

/*
 * ===========================================================================
 * PRODUCTION $Log: density_map.hpp,v $
 * PRODUCTION Revision 1000.0  2004/06/01 19:54:41  gouriano
 * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.1
 * PRODUCTION
 * ===========================================================================
 */

#ifndef GUI_UTILS___DENSITY_GRAPH__HPP
#define GUI_UTILS___DENSITY_GRAPH__HPP

/*  $Id: density_map.hpp,v 1000.0 2004/06/01 19:54:41 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:  Mike DiCuccio
 *
 * File Description:
 *
 */

#include <corelib/ncbistd.hpp>
#include <objmgr/bioseq_handle.hpp>
#include <objmgr/feat_ci.hpp>
#include <objmgr/align_ci.hpp>
#include <objmgr/annot_ci.hpp>
#include <objtools/alnmgr/alnmap.hpp>
#include <objtools/alnmgr/alnmix.hpp>

#include <util/range.hpp>
#include <deque>


/** @addtogroup GUI_UTILS
 *
 * @{
 */


BEGIN_NCBI_SCOPE


/// Run iterator.  iterate through runs of equal values in the bins.


template <typename CntType>
class CDenMapRunIterator {
public:
    typedef vector<CntType> container_type;
    typedef typename container_type::size_type   position_type;
        
    CDenMapRunIterator(position_type here, const container_type& bins, TSeqPos start, TSeqPos window)
    : m_Bins(bins), m_Pos(here), m_RunLength(x_CalcRunLength()), m_Start(start), m_Window(window)
    {}

    CntType       Advance();
        
    position_type GetPosition() const { return m_Pos; }
    position_type GetRunLength() const{ return m_RunLength; }
    TSeqPos       GetSeqPosition() const;
    TSeqPos       GetSeqRunLength() const;
    
    CntType       GetValue() const { return Valid() ? m_Bins[m_Pos]: 0 ; }
    
    bool operator ==(CDenMapRunIterator<CntType>& rhs) const { return m_Pos == rhs.m_Pos; }
    
    operator bool() const { return Valid(); }
    bool        Valid() const { return  m_Pos >= 0 &&  m_Pos < m_Bins.size() && m_RunLength != 0; }
private:

    
    position_type   x_CalcRunLength() const;
    
    const container_type&   m_Bins;
    position_type       m_Pos;
    position_type       m_RunLength;
    
    const TSeqPos       m_Start;
    const TSeqPos       m_Window;
    
};

template <typename CntType>
TSeqPos CDenMapRunIterator<CntType>::GetSeqPosition() const
{
    return GetPosition() * m_Window + m_Start; 
}

template <typename CntType>
TSeqPos CDenMapRunIterator<CntType>::GetSeqRunLength() const 
{ 
    return GetRunLength() * m_Window; 
}


template <typename CntType>
typename CDenMapRunIterator<CntType>::position_type 
CDenMapRunIterator<CntType>::x_CalcRunLength() const
{
    if (m_Pos >= 0 &&  m_Pos < m_Bins.size()) {
        CntType this_value = m_Bins[m_Pos];
        position_type i;
        for (i = m_Pos + 1; 
             i < m_Bins.size()  &&  m_Bins[i] == this_value;
             ++i) {}

        return i - m_Pos;
    }
    return 0;
}

template <typename CntType>
CntType
CDenMapRunIterator<CntType>::Advance() {
    m_Pos += m_RunLength;
    m_RunLength = x_CalcRunLength();
    return GetValue();
}


///
/// class CDensityMap generates a low-resolution view of a set of features.
/// The features to be processed are defined by the SAnnotSelector object.
/// The return value is a vector of counts, one for each of the bins
/// defined by the window width and spanning the range from start to stop.
/// If start and stop are both zero, then the entire sequence is evaluated.
///


template <typename CntType, typename Accum  = plus<CntType> >
class CDensityMap {
public:
    typedef vector<CntType> container_type;
    typedef typename container_type::const_iterator const_iterator;
    typedef typename container_type::iterator       iterator;
    typedef CDenMapRunIterator<CntType>             runlen_iterator;
    
    CDensityMap(TSeqPos start, TSeqPos stop, TSeqPos window = 1,
                Accum ac = Accum());
    CDensityMap(const objects::CBioseq_Handle& handle, TSeqPos window = 1,
                Accum ac = Accum());
    
    void AddRange(TSeqRange range, CntType score = 1, bool expand = false);
    
    /// All features on this bioseq selected by sel in the range of this.
    CntType AddFeatures(const objects::CBioseq_Handle& handle,
                        objects::SAnnotSelector sel);
    
    /// All alignments on this bioseq selected by sel in the range of this.
    CntType AddAlignments(const objects::CBioseq_Handle& handle, 
                          objects::SAnnotSelector sel);

    /// All alignments in a given annotation on this bioseq within the range of this.
    CntType AddAlignments(const objects::CBioseq_Handle& handle, 
                          const objects::CSeq_annot& seq_annot);

    void 		Clear() { m_Max = 0; fill(m_Bins.begin(), m_Bins.end(), 0); }
    TSeqPos     GetStart() const { return m_Range.GetFrom(); }
    TSeqPos     GetStop() const { return m_Range.GetTo(); }
    TSeqRange   GetRange() const { return m_Range; }
    TSeqPos     GetWindow() const { return m_Window; }
    size_t      GetBins() const { return m_Bins.size(); }
    CntType     GetMax() const { return m_Max; }
    Accum       GetAccum() const { return m_AccumFunc; }
    
    /// extend our density map to cover the sequence position stop.
    /// can only be used to extend to the right, that is only Stop is affected, not Start.
    void            ExtendTo(TSeqPos stop);
    
    const_iterator  begin() const { return m_Bins.begin(); }
    const_iterator  end() const { return m_Bins.end(); }
    iterator        begin()  { return m_Bins.begin(); }
    iterator        end()  { return m_Bins.end(); }
    
    runlen_iterator RunLenBegin() const { return RunLenIterator(0); }
    runlen_iterator RunLenIterator(typename container_type::size_type n) const 
        { return runlen_iterator(n, m_Bins, m_Range.GetFrom(), m_Window); }
    
    CntType operator[](typename container_type::size_type n) { return m_Bins[n]; }

    /// OLD static method. Use AddFeatures method instead.
    /// retrieve a density map.  The return value is the maximum value
    /// in the density graph.
    static TSeqPos GetDensityMap(const objects::CBioseq_Handle& handle,
                                 TSeqPos start, TSeqPos stop,
                                 TSeqPos window,
                                 objects::SAnnotSelector sel,
                                 vector<TSeqPos>& density);

private:
    /// given the range and window size, how many bins should there be?
    size_t  x_CalcNbins() { return (m_Range.GetTo() - m_Range.GetFrom() + m_Window) / m_Window; }
    /// convert from sequence coords to a bin number.
    size_t  x_BinN(TSeqPos p) { return (p - m_Range.GetFrom())/m_Window; }

    /// closed range on a sequence this covers.
    TSeqRange       m_Range;

    /// coordinates per bin.
    TSeqPos         m_Window;

    /// maximum Count accumulated in the bins so far.
    CntType         m_Max;
    
    /// Where we actually keep the accumulated counts/scores/whatever.
    container_type  m_Bins;
    
    Accum           m_AccumFunc;
};


template <typename CntType, typename Accum >
CDensityMap<CntType, Accum>::CDensityMap(TSeqPos start, 
                                         TSeqPos stop,
                                         TSeqPos window /* = 1 */, 
                                         Accum ac /*= Accum() */)
:   m_Range(start, stop), 
    m_Window(window < 1 ? 1 : window), 
    m_Max(0), 
    m_Bins(x_CalcNbins()), 
    m_AccumFunc(ac)
{
}


template <typename CntType, typename Accum >
CDensityMap<CntType, Accum>::CDensityMap(const objects::CBioseq_Handle& handle, 
                                         TSeqPos window /* = 1 */, 
                                         Accum ac /*= Accum() */)
:   m_Range(0, handle.GetSeqVector().size()), 
    m_Window(window < 1 ? 1 : window),  
    m_Max(0), 
    m_Bins(x_CalcNbins()), 
    m_AccumFunc(ac)
{
}


template <typename CntType, typename Accum >
void  CDensityMap<CntType, Accum>::ExtendTo(TSeqPos stop) 
{
    if (stop > GetStop()) {
        m_Range.SetTo(stop);
        m_Bins.resize(x_CalcNbins(), 0);
    }
}


template <typename CntType, typename Accum >
void CDensityMap<CntType, Accum>::AddRange(TSeqRange range, CntType score, bool expand)
{   
    //_ASSERT(range.GetFrom() <= range.GetTo());
    if (range.GetFrom() > range.GetTo()) {
        range = TSeqRange(range.GetTo(), range.GetFrom());
    }
    
    if ( expand ) {
        ExtendTo( range.GetTo());
    }
    
    TSeqRange usable_range(m_Range.IntersectionWith(range));
    if (usable_range.Empty()) {
        return;
    }
    
    size_t begin_bin = x_BinN(usable_range.GetFrom());
    size_t end_bin   = x_BinN(usable_range.GetTo()) + 1;
    for (size_t i = begin_bin;  i < end_bin ;  ++i) {
        CntType new_val = m_AccumFunc(m_Bins[i], score);
        m_Bins[i] = new_val;
        m_Max = max(m_Max, new_val);
    }
}


template <typename CntType, typename Accum >
CntType CDensityMap<CntType, Accum>::AddFeatures(
    const objects::CBioseq_Handle& handle, 
    objects::SAnnotSelector sel)
{
    TSeqPos start(GetStart());
    TSeqPos stop(GetStop());

   // grab a feature iterator for our range
    objects::CFeat_CI feat_iter(handle, start, stop, sel);
    if (feat_iter.GetSize() == 0) {
        return 0;
    }

    for (;  feat_iter;  ++feat_iter) {
        objects::CMappedFeat feat = *feat_iter;
        TSeqRange range = feat.GetLocation().GetTotalRange();
        AddRange(range, 1, false);
    }

    return GetMax();

}
    

template <typename CntType, typename Accum >
CntType CDensityMap<CntType, Accum>::AddAlignments(
    const objects::CBioseq_Handle& handle, 
    objects::SAnnotSelector sel)
{
    TSeqPos start(GetStart());
    TSeqPos stop(GetStop());