splign_hit.cpp

ncbi源码

}


bool CHit::IsConsistent() const
{
    return m_ai[0] <= m_ai[1] && m_ai[2] <= m_ai[3];
}


void CHit::SetEnvelop()
{
    m_ai[0] = min(m_an[0],m_an[1]);    m_ai[1] = max(m_an[0],m_an[1]);
    m_ai[2] = min(m_an[2],m_an[3]);    m_ai[3] = max(m_an[2],m_an[3]);
}


void CHit::Move(unsigned char i, int pos_new, bool prot2nucl)
{
    bool strand_plus = IsStraight();
    // for minus strand, flip the hit if necessary
    if(!strand_plus && m_an[0] > m_an[1]) {
        int n = m_an[0]; m_an[0] = m_an[1]; m_an[1] = n;
        n = m_an[2]; m_an[2] = m_an[3]; m_an[3] = n;
    }
    
    // i --> m_ai, n --> m_an
    unsigned char n = 0;
    switch(i) {
        case 0: n = 0; break; 
        case 1: n = 1; break; 
        case 2: n = strand_plus? 2: 3; break; 
        case 3: n = strand_plus? 3: 2; break; 
    }
    
    int shift = pos_new - m_an[n], shiftQ, shiftS;
    if(shift == 0) return;

    // figure out shifts on query and subject
    if(prot2nucl) {
        if(n < 2) {
            shiftQ = shift;
            shiftS = 3*shiftQ;
        }
        else {
            if(shift % 3) {
                shiftQ = (shift > 0)? (shift/3 + 1): (shift/3 - 1);
                shiftS = 3*shiftQ;
            }
            else {
                shiftQ = shift / 3;
                shiftS = shift;
            }
        }
    }
    else {
        shiftQ = shiftS = shift;
    }

    // make sure hit ratio doesn't change (todo: prot2nucl)
    if(!prot2nucl) {
        bool plus = shift >= 0;
        unsigned shift_abs = plus? shift: -shift;
        if(n < 2) {
            double D = double(shift_abs) / (m_ai[1] - m_ai[0] + 1);
            shiftQ = shift;
            shiftS = int(ceil(D*(m_ai[3] - m_ai[2] + 1)));
            if(!plus) shiftS = -shiftS;
        }
        else {
            double D = double(shift_abs) / (m_ai[3] - m_ai[2] + 1);
            shiftS = shift;
            shiftQ = int(ceil(D*(m_ai[1] - m_ai[0] + 1)));
            if(!plus) shiftQ = -shiftQ;
        }
    }


    if(strand_plus) {
        m_an[n % 2] += shiftQ; m_an[n % 2 + 2] += shiftS;
        m_ai[i % 2] += shiftQ; m_ai[i % 2 + 2] += shiftS;
    }
    else {
        switch(i) {
            case 0: m_an[0] = m_ai[0] += shiftQ;
                    m_an[2] = m_ai[3] -= shiftS;
                    break;
            case 1: m_an[1] = m_ai[1] += shiftQ;
                    m_an[3] = m_ai[2] -= shiftS;
                    break;
            case 2: m_an[3] = m_ai[2] += shiftS;
                    m_an[1] = m_ai[1] -= shiftQ;
                    break;
            case 3: m_an[2] = m_ai[3] += shiftS;
                    m_an[0] = m_ai[0] -= shiftQ;
                    break;
        }
    }
}


bool CHit::Inside(const CHit& b) const
{
    bool b0 = b.m_ai[0] <= m_ai[0] && m_ai[1] <= b.m_ai[1];
    if(!b0) return false;
    bool b1 = b.m_ai[2] <= m_ai[2] && m_ai[3] <= b.m_ai[3];
    return b1;
}

// This function must be called once after a series of Move() calls.
// Changes are made based on previous values of the parameters.
void CHit::UpdateAttributes()
{
    int nLength0 = m_Length;
    m_anLength[0] = m_ai[1] - m_ai[0] + 1;
    m_anLength[1] = m_ai[3] - m_ai[2] + 1;
    m_Length = max(m_anLength[0], m_anLength[1]);
    double dk = nLength0? double(m_Length) / nLength0: 1.0;
    m_Score *= dk;
    m_Gaps = int(m_Gaps*dk);
    m_MM = int(m_MM*dk);
    // do not change m_Value, m_Idnty
}

ostream& operator << (ostream& os, const CHit& hit) {
    const int nLen = min(hit.m_anLength[0], hit.m_anLength[1]);
    return os
       << hit.m_Query << '\t' << hit.m_Subj << '\t'
       << hit.m_Idnty   << '\t' << nLen          << '\t'
       << hit.m_MM      << '\t' << hit.m_Gaps   << '\t'
       << hit.m_an[0]    << '\t' << hit.m_an[1]   << '\t'
       << hit.m_an[2]    << '\t' << hit.m_an[3]   << '\t'
       << hit.m_Value   << '\t' << hit.m_Score;
}


bool hit_same_order::operator() (const CHit& hm, const CHit& h0) const
{
    bool bStraightM = hm.IsStraight();
    bool bStraight0 = h0.IsStraight();
    bool b0 = bStraightM && bStraight0;
    bool b1 = !bStraightM && !bStraight0;
    if(b0 || b1)
    {
        int vm [4]; copy(hm.m_an, hm.m_an + 4, vm);
        int v0 [4]; copy(h0.m_an, h0.m_an + 4, v0);
        if(b1)
        {
                // we need to change direction for one axis
            double dmin1 = 1e38, dmax1 = -dmin1;
            int i = 0;
            for(i=2; i<4; i++)
            {
                if(dmin1 > vm[i]) dmin1 = vm[i];
                   if(dmin1 > v0[i]) dmin1 = v0[i];
                if(dmax1 < vm[i]) dmax1 = vm[i];
                   if(dmax1 < v0[i]) dmax1 = v0[i];
            }
            for(i=2; i<4; i++)
            {
                vm[i] = int(-vm[i] + dmin1 + dmax1);
                v0[i] = int(-v0[i] + dmin1 + dmax1);
            }
                // now u-turn hits if necessary
            if(vm[0] > vm[1] && vm[2] > vm[3])
            {
                int n = vm[0]; vm[0] = vm[1]; vm[1] = n;
                n = vm[2]; vm[2] = vm[3]; vm[3] = n;
            }
            if(v0[0] > v0[1] && v0[2] > v0[3])
            {
                int n = v0[0]; v0[0] = v0[1]; v0[1] = n;
                n = v0[2]; v0[2] = v0[3]; v0[3] = n;
            }
        }
        double dx = .1*0.5*(abs(vm[1] - vm[0]) + abs(v0[1] - v0[0]));
        if(vm[0] > v0[1] - dx) // left bottom
            return vm[2] > v0[3] - dx? true: false;
        else
        if(vm[2] > v0[3] - dx) // left top
            return vm[0] > v0[1] - dx? true: false;
        else
        if(vm[1] < v0[0] + dx) // right bottom
            return vm[3] < v0[2] + dx? true: false;
        else
        if(vm[3] < v0[2] + dx) // right top
            return vm[1] < v0[0] + dx? true: false;
        else return false;
    }
    else
        return false;   // different strands and same order are incompatible
}

bool hit_not_same_order::operator()(const CHit& hm, const CHit& h0) const
{
    return !hit_same_order::operator() (hm,h0);
}

// Calculates distance btw. two hits on query or subject
int CHit::GetHitDistance(const CHit& h1, const CHit& h2, int nWhere)
{
    int i0 = -1, i1 = -1;
    switch(nWhere)
    {
        case 0: // query
            i0 = 0; i1 = 1;
            break;
        case 1: // subject
            i0 = 2; i1 = 3;
            break;
        default:
            return -2; // wrong parameter
    }
    if(h1.m_ai[i1] < h2.m_ai[i0])  // no intersection
        return h2.m_ai[i0] - h1.m_ai[i1];
    else
    if(h2.m_ai[i1] < h1.m_ai[i0])  // no intersection
        return h1.m_ai[i0] - h2.m_ai[i1];
    else // overlapping
        return 0;
}

END_NCBI_SCOPE


/*
* ===========================================================================
*
* $Log: splign_hit.cpp,v $
* Revision 1000.0  2004/06/01 18:12:48  gouriano
* PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.8
*
* Revision 1.8  2004/05/24 16:13:57  gorelenk
* Added PCH ncbi_pch.hpp
*
* Revision 1.7  2004/05/10 16:33:52  kapustin
* Report unexpected strands in CHit(CSeq_align&). Calculate length and identity based on alignment data.
*
* Revision 1.6  2004/05/06 17:43:27  johnson
* CHit(CSeq_align&) now indicates minus strand via swapping subject
* coordinates (*not* query coordinates)
*
* Revision 1.5  2004/05/04 15:23:45  ucko
* Split splign code out of xalgoalign into new xalgosplign.
*
* Revision 1.4  2004/05/03 15:23:08  johnson
* Added CHit constructor that takes a CSeq_align
*
* Revision 1.3  2004/04/23 14:37:44  kapustin
* *** empty log message ***
*
* 
* ===========================================================================
*/