parse.cpp

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: parse.cpp,v $
 * PRODUCTION Revision 1000.2  2004/06/01 18:08:37  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.3
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: parse.cpp,v 1000.2 2004/06/01 18:08:37 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:  Alexandre Souvorov
 *
 * File Description:
 *
 */

#include <ncbi_pch.hpp>
#include "gene_finder.hpp"

BEGIN_NCBI_SCOPE

template<class T> inline void PushState(ParseVec<T>* vecp, int strand, int point)
{
    for (int phase = 0;  phase < 3;  ++phase) {
        vecp[phase].push_back(T(strand,phase,point));
        ++vecp[phase].num;
    }
}


template<class T> inline void PushState(ParseVec<T>& vec, int strand, int point)
{
    vec.push_back(T(strand,point));
    ++vec.num;
}


inline void PushState(ParseVec<LastExon>& vec, int strand, int point)
{
    vec.push_back(LastExon(strand,2,point));
    ++vec.num;
}


inline void PushState(ParseVec<FirstExon>& vec, int strand, int point)
{
    vec.push_back(FirstExon(strand,0,point));
    ++vec.num;
}


template<class L, class R> inline 
bool TooFar(const L& left, const R& right, double score) 
{ 
    if (left.MScore() == BadScore) {
        return true;
    }
    int len = right.Stop() - left.Stop();
    return (len > TooFarLen  &&  score + left.MScore() < right.Score());
}


struct RState
{
    RState(const HMM_State& l, const HMM_State& r)
    {
        lsp = r.LeftState();
        rsp = const_cast<HMM_State*>(&r);
        rsp->UpdateLeftState(l);
    }
    ~RState() { rsp->UpdateLeftState(*lsp); } 

    const HMM_State* lsp;
    HMM_State* rsp;
};

template<class Left, class Right> inline
bool EvaluateNewScore(const Left& left, const Right& right, double& rscore, 
                      bool& openrgn, bool& inalign)
{

    rscore = BadScore;

    RState rst(left,right);

    int len = right.Stop() - left.Stop();
    if (len > right.MaxLen()) {
        return false;
    }
    if ( !right.NoRightEnd()  &&  len < right.MinLen() ) {
        return true;
    }

    double scr, score = 0;
    if (left.isPlus()) {
        scr = left.BranchScore(right);
        if (scr == BadScore) {
            return true;
        }
    } else {
        scr = right.BranchScore(left);
        if (scr == BadScore) {
            return true;
        }
        scr += right.DenScore() - left.DenScore();
    }
    score += scr;

    // this check is frame - dependent and MUST be after BranchScore
    if (right.StopInside()) {
        return false;
    }     

    if (right.NoRightEnd()) {
        scr = right.ClosingLengthScore();
    } else {
        scr = right.LengthScore();
    }
    if (scr == BadScore) {
        return true;
    }
    score += scr;

    scr = right.RgnScore();
    if (scr == BadScore) {
        return true;
    }
    score += scr;

    if ( !right.NoRightEnd() ) {
        scr = right.TermScore();
        if (scr == BadScore) {
            return true;
        }
        score += scr;
    }

    openrgn = right.OpenRgn();
    inalign = right.InAlignment();
    rscore = score;
    return true;
}


class RightIsExon {};
class RightIsNotExon {};

template<class L, class R> // template for all right exons
inline bool ForwardStep(const L& left, R& right, RightIsExon rie)
{
    double score;
    bool openrgn, inalign;
    if ( !EvaluateNewScore(left,right,score,openrgn,inalign) ) {
        return false;
    } else if (score == BadScore) {
        return true;
    }

    if (left.Score() != BadScore  &&  openrgn  &&  !inalign) {
        double scr = score + left.Score();
        if (scr > right.Score()) {
            right.UpdateLeftState(left);
            right.UpdateScore(scr);
        }
    }

    if ( !openrgn ) {
        return false;
    }

    return true;
}


template<class L, class R> // template for all right non - exons
inline bool ForwardStep(const L& left, R& right, RightIsNotExon rine)
{
    double score;
    bool openrgn, inalign;
    if ( !EvaluateNewScore(left,right,score,openrgn,inalign) ) {
        return false;
    } else if (score == BadScore) {
        return true;
    }

    if (left.Score() != BadScore  &&  openrgn  &&  !inalign) {
        double scr = score + left.Score();
        if (scr > right.Score()) {
            right.UpdateLeftState(left);
            right.UpdateScore(scr);
        }
    }

    if ( !openrgn  ||  TooFar(left, right, score) ) {
        return false;
    }

    return true;
}


template<class L, class R> // template for all right exons
inline void MakeStep(ParseVec<L>& lvec, ParseVec<R>& rvec, RightIsExon rie)
{
    if (lvec.num < 0) {
        return;
    }
    int i = lvec.num;
    if (lvec[i].Stop() == rvec[rvec.num].Stop()) {
        --i;
    }
    while (i >= 0  &&  ForwardStep(lvec[i--],rvec[rvec.num],rie)) {
    }

    if (rvec.num > 0) {
        rvec[rvec.num].UpdatePrevExon(rvec[rvec.num - 1]);
    }
}


template<class L, class R> // template for all right non - exons
inline void MakeStep(ParseVec<L>& lvec, ParseVec<R>& rvec, RightIsNotExon rine)
{
    if (lvec.num < 0) {
        return;
    }
    R& right = rvec[rvec.num];
    int i = lvec.num;
    int rlimit = right.Stop();
    if (lvec[i].Stop() == rlimit) {
        --i;
    }
    int nearlimit = max(0,rlimit - TooFarLen);
    while (i >= 0  &&  lvec[i].Stop() >= nearlimit) {
        if ( !ForwardStep(lvec[i--],right,rine) ) {
            return;
        }
    }
    if (i < 0) {
        return;
    }
    for (const L* p = &lvec[i];  p !=0;  p = p->PrevExon()) {
        if ( !ForwardStep(*p,right,rine) ) {
            return;
        }
    }
}

// R / L indicates that left / right parameter is a 3 - element array

template<class L, class R, class D> 
inline void MakeStepR(L& lvec, R rvec, D d)
{
    for (int phr = 0;  phr < 3;  ++phr) {
        MakeStep(lvec, rvec[phr], d);
    }
}


template<class L, class R, class D> 
inline void MakeStepL(L lvec, R& rvec, D d)
{
    for (int phl = 0;  phl < 3;  ++phl) {
        MakeStep(lvec[phl], rvec, d);
    }
}


template<class L, class R, class D> 
inline void MakeStepLR(L lvec, R rvec, D d)
{
    for (int phl = 0;  phl < 3;  ++phl) {
        for (int phr = 0;  phr < 3;  ++phr)
        {
            MakeStep(lvec[phl], rvec[phr], d);
        }
    }
}


template<class L, class R, class D> 
inline void MakeStepLR(L lvec, R rvec, D d, int shift)
{
    for (int phl = 0;  phl < 3;  ++phl) {
        int phr = (shift + phl)%3;
        MakeStep(lvec[phl], rvec[phr], d);
    }
}


Parse::Parse(const SeqScores& ss) : seqscr(ss)
{
    try {
        igplus.reserve (seqscr.StartNumber(Plus) + 1);
        igminus.reserve(seqscr.StopNumber (Minus) + 1);
        feminus.reserve(seqscr.StartNumber(Minus) + 1);
        leplus.reserve (seqscr.StopNumber (Plus) + 1);
        seplus.reserve (seqscr.StopNumber (Plus) + 1);
        seminus.reserve(seqscr.StartNumber(Minus) + 1);
        for (int phase = 0;  phase < 3;  ++phase) {
            inplus[phase ].reserve(seqscr.AcceptorNumber(Plus) + 1);
            inminus[phase].reserve(seqscr.DonorNumber   (Minus) + 1);
            feplus[phase ].reserve(seqscr.DonorNumber   (Plus) + 1);
            ieplus[phase ].reserve(seqscr.DonorNumber   (Plus) + 1);
            ieminus[phase].reserve(seqscr.AcceptorNumber(Minus) + 1);
            leminus[phase].reserve(seqscr.AcceptorNumber(Minus) + 1);
        }
    }
    catch(bad_alloc&) {
        NCBI_THROW(CGnomonException, eGenericError,
                   "Not enough memory in Parse");
    }

    int len = seqscr.SeqLen();

    RightIsExon rie;
    RightIsNotExon rine;

    for (int i = 0;  i < len;  ++i) {
        if (seqscr.AcceptorScore(i,Plus) != BadScore) {
            PushState (inplus,Plus,i);
            MakeStepLR(ieplus,inplus,rine,1);
            MakeStepLR(feplus,inplus,rine,1);
        }

        if (seqscr.AcceptorScore(i,Minus) != BadScore) {
            PushState (ieminus,Minus,i);
            PushState (leminus,Minus,i);
            MakeStepLR(inminus,ieminus,rie);
            MakeStepR(igminus,leminus,rie);
        }

        if (seqscr.DonorScore(i,Plus) != BadScore) {
            PushState (feplus,Plus,i);
            PushState (ieplus,Plus,i);
            MakeStepLR(inplus,ieplus,rie);
            MakeStepR(igplus,feplus,rie);
        }

        if (seqscr.DonorScore(i,Minus) != BadScore) {
            PushState (inminus,Minus,i);
            MakeStepLR(leminus,inminus,rine,0);
            MakeStepLR(ieminus,inminus,rine,0);
        }

        if (seqscr.StartScore(i,Plus) != BadScore) {
            PushState(igplus,Plus,i);
            MakeStep (seplus,igplus,rine);
            MakeStep (seminus,igplus,rine);
            MakeStep (leplus,igplus,rine);
            MakeStep (feminus,igplus,rine);
        }

        if (seqscr.StartScore(i,Minus) != BadScore) {
            PushState(feminus,Minus,i);
            PushState(seminus,Minus,i);
            MakeStepL(inminus,feminus,rie);
            MakeStep (igminus,seminus,rie);
        }

        if (seqscr.StopScore(i,Plus) != BadScore) {
            PushState(leplus,Plus,i);
            PushState(seplus,Plus,i);
            MakeStepL(inplus,leplus,rie);
            MakeStep (igplus,seplus,rie);
        }

        if (seqscr.StopScore(i,Minus) != BadScore) {
            PushState(igminus,Minus,i);
            MakeStep (seplus,igminus,rine);
            MakeStep (seminus,igminus,rine);
            MakeStep (leplus,igminus,rine);
            MakeStep (feminus,igminus,rine);
        }
    }

    PushState (inplus,Plus,-1);
    MakeStepLR(ieplus,inplus,rine,1);
    MakeStepLR(feplus,inplus,rine,1);

    PushState (ieminus,Minus,-1);
    MakeStepLR(inminus,ieminus,rie);

    PushState(leminus,Minus,-1);
    MakeStepR(igminus,leminus,rie);

    PushState (ieplus,Plus,-1);
    MakeStepLR(inplus,ieplus,rie);

    PushState(feplus,Plus,-1);
    MakeStepR(igplus,feplus,rie);

    PushState (inminus,Minus,-1);
    MakeStepLR(leminus,inminus,rine,0);
    MakeStepLR(ieminus,inminus,rine,0);

    PushState(igplus,Plus,-1);
    MakeStep (seplus,igplus,rine);
    MakeStep (seminus,igplus,rine);
    MakeStep (leplus,igplus,rine);
    MakeStep (feminus,igplus,rine);

    PushState(feminus,Minus,-1);
    MakeStepL(inminus,feminus,rie);

    PushState(seminus,Minus,-1);
    MakeStep (igminus,seminus,rie);

    PushState(leplus,Plus,-1);
    MakeStepL(inplus,leplus,rie);

    PushState(seplus,Plus,-1);
    MakeStep(igplus,seplus,rie);

    PushState(igminus,Minus,-1);
    MakeStep (seplus,igminus,rine);
    MakeStep (seminus,igminus,rine);
    MakeStep (leplus,igminus,rine);
    MakeStep (feminus,igminus,rine);

    igplus.back().UpdateScore(AddProbabilities(igplus.back().Score(),igminus.back().Score()));

    const HMM_State* p = &igplus.back();
    if (feminus.back().Score() > p->Score()) {
        p = &feminus.back();
    }
    if (leplus.back().Score() > p->Score()) {
        p = &leplus.back();
    }
    if (seplus.back().Score() > p->Score()) {
        p = &seplus.back();
    }
    if (seminus.back().Score() > p->Score()) {
        p = &seminus.back();
    }
    for (int ph = 0;  ph < 3;  ++ph) {
        if (inplus[ph].back().Score() > p->Score()) {
            p = &inplus[ph].back();
        }
        if (inminus[ph].back().Score() > p->Score()) {
            p = &inminus[ph].back();
        }
        if (ieplus[ph].back().Score() > p->Score()) {
            p = &ieplus[ph].back();
        }
        if (ieminus[ph].back().Score() > p->Score()) {
            p = &ieminus[ph].back();