models.cpp

ncbi源码

    }
    for (int i = 1;  i < len;  ++i) {
        notining[i] = max(notining[i - 1],notining[i]);
    }

    const int stpT = 1, stpTA = 2, stpTG = 4;

    for (int strand = 0;  strand < 2;  ++strand) {
        IVec& s = seq[strand];

        anum[strand] = 0;
        dnum[strand] = 0;
        sttnum[strand] = 0;
        stpnum[strand] = 0;

        if (strand == Plus) {
            for (int i = 0;  i < len;  ++i) {
                ascr[strand][i] = max(ascr[strand][i],acceptor.Score(s,i));
                if (ascr[strand][i] != BadScore) {
                    if (s[i + 1] == nA  &&  s[i + 2] == nA) {
                        asplit[strand][0][i] |= stpT;
                    }
                    if (s[i + 1] == nA  &&  s[i + 2] == nG) {
                        asplit[strand][0][i] |= stpT;
                    } 
                    if (s[i + 1] == nG  &&  s[i + 2] == nA) {
                        asplit[strand][0][i] |= stpT;
                    }
                    if (s[i + 1] == nA) {
                        asplit[strand][1][i] |= stpTA|stpTG;
                    }
                    if (s[i + 1] == nG) {
                        asplit[strand][1][i] |= stpTA;
                    }
                }
                dscr[strand][i] = max(dscr[strand][i],donor.Score(s,i));
                if (dscr[strand][i] != BadScore) {
                    if (s[i] == nT) {
                        dsplit[strand][0][i] |= stpT;
                    }
                    if (s[i - 1] == nT  &&  s[i] == nA) {
                        dsplit[strand][1][i] |= stpTA;
                    }
                    if (s[i - 1] == nT  &&  s[i] == nG) {
                        dsplit[strand][1][i] |= stpTG;
                    }
                }
                sttscr[strand][i] = start.Score(s,i);
                stpscr[strand][i] = stop.Score(s,i);
                if (ascr[strand][i] != BadScore) {
                    ++anum[strand];
                }
                if (dscr[strand][i] != BadScore) {
                    ++dnum[strand];
                }
                if (sttscr[strand][i] != BadScore) {
                    ++sttnum[strand];
                }
            }
        } else {
            for (int i = 0;  i < len;  ++i) {
                int ii = len - 2-i;   // extra -1 because ii is point on the "right"
                ascr[strand][i] = max(ascr[strand][i],acceptor.Score(s,ii));
                if (ascr[strand][i] != BadScore) {
                    if (s[ii + 1] == nA  &&  s[ii + 2] == nA) {
                        asplit[strand][0][i] |= stpT;
                    }
                    if (s[ii + 1] == nA  &&  s[ii + 2] == nG) {
                        asplit[strand][0][i] |= stpT;
                    } 
                    if (s[ii + 1] == nG  &&  s[ii + 2] == nA) {
                        asplit[strand][0][i] |= stpT;
                    }
                    if (s[ii + 1] == nA) {
                        asplit[strand][1][i] |= stpTA|stpTG;
                    }
                    if (s[ii + 1] == nG) {
                        asplit[strand][1][i] |= stpTA;
                    }
                }
                dscr[strand][i] = max(dscr[strand][i],donor.Score(s,ii));
                if (dscr[strand][i] != BadScore) {
                    if (s[ii] == nT) {
                        dsplit[strand][0][i] |= stpT;
                    }
                    if (s[ii - 1] == nT  &&  s[ii] == nA) {
                        dsplit[strand][1][i] |= stpTA;
                    }
                    if (s[ii - 1] == nT  &&  s[ii] == nG) {
                        dsplit[strand][1][i] |= stpTG;
                    }
                }
                sttscr[strand][i] = start.Score(s,ii);
                stpscr[strand][i] = stop.Score(s,ii);
                if (ascr[strand][i] != BadScore) {
                    ++anum[strand];
                }
                if (dscr[strand][i] != BadScore) {
                    ++dnum[strand];
                }
            }
        }
    }       

    for (CClusterSet::ConstIt it_cls = cls_local.begin();  it_cls != cls_local.end();  ++it_cls) {
        for (CCluster::ConstIt it = it_cls->begin();  it != it_cls->end();  ++it) {
            const AlignVec& algn = *it;
            int strand = algn.Strand();
            int type = algn.Type();

            if (type != AlignVec::Prot) {
                IVec& s = seq[Plus];
                IVec mRNA;
                for (int k = 0;  k < algn.size();  ++k) {
                    int exona = algn[k].first;
                    int exonb = algn[k].second;
                    copy(&s[exona],&s[exonb]+1,back_inserter(mRNA));
                }

                if (strand == Plus) {
                    int shft = algn.front().first;
                    for (int k = 0;  k < algn.size() - 1;  ++k) {
                        int exonb = algn[k].second;

                        for (int i = max(0,exonb - 2);  i <= exonb;  ++i)    // if i == exonb, the stop is in the next "exon"
                        {
                            stpscr[strand][i] = stop.Score(mRNA,i - shft);
                        }

                        int exona = algn[k + 1].first;   //next exon
                        shft += exona - exonb - 1;         //intron length

                        for (int i = exona - 1;  i <= min(len - 1,exona + 1);  ++i) {
                            sttscr[strand][i] = start.Score(mRNA,i - shft);
                        }
                    }
                } else {
                    reverse(mRNA.begin(),mRNA.end());
                    for (int i = 0;  i < mRNA.size();  ++i) mRNA[i] = toMinus[mRNA[i]];

                    int shft = algn.back().second - 1;
                    for (int k = algn.size() - 1;  k > 0;  --k) {
                        int exona = algn[k].first;

                        for (int i = exona - 1;  i <= min(len - 1,exona + 1);  ++i) {
                            stpscr[strand][i] = stop.Score(mRNA,shft - i);
                        }

                        int exonb = algn[k - 1].second;  //prev exon
                        shft -= exona - exonb - 1;         //intron length

                        for (int i = max(0,exonb - 2);  i <= exonb;  ++i) {
                            sttscr[strand][i] = start.Score(mRNA,shft - i);
                        }
                    }
                }
            }
        }
    }

    for (CClusterSet::ConstIt it_cls = cls_local.begin();  it_cls != cls_local.end();  ++it_cls) {
        for (CCluster::ConstIt it = it_cls->begin();  it != it_cls->end();  ++it) {
            const AlignVec& algn = *it;
            int strand = algn.Strand();
            int type = algn.Type();

            if (type == AlignVec::Prot) {
                for (int k = 0;  k < algn.size();  ++k)   // suppressing STOPS in prot align
                {
                    int exona = algn[k].first;
                    int exonb = algn[k].second;
                    for (int i = exona;  i < exonb;  i += 3) {
                        if (strand == Plus) {
                            if (i > 0) {
                                stpscr[strand][i - 1] = BadScore;
                            }  // score on last coding base
                        } else {
                            if (i + 2 < len - 1) {
                                stpscr[strand][i + 2] = BadScore;
                            }  // score on last intergenic(first stop) base
                        }
                    }
                }
            }
        }
    }

    for (int strand = 0;  strand < 2;  ++strand) {
        for (int i = 0;  i < len;  ++i) {
            if (sttscr[strand][i] != BadScore) {
                ++sttnum[strand];
            }

            if (stpscr[strand][i] != BadScore) {
                if (strand == Plus)
                {
                    int& lstp = laststop[strand][2 - i%3][i + 3];
                    lstp = max(i + 1,lstp);
                    ++stpnum[strand];
                } else {
                    int& lstp = laststop[strand][i%3][i];
                    lstp = max(i - 2,lstp);
                    ++stpnum[strand];
                }
            }
        }
    }


    for (int strand = 0;  strand < 2;  ++strand) {
        IVec& s = seq[strand];

        for (int i = 0;  i < len;  ++i) {
            int ii = strand == Plus ? i : len - 1-i;

            double score = ncdr.Score(s,ii);
            if (score == BadScore) {
                score = 0;
            }
            ncdrscr[strand][i] = score;
            if (i > 0) {
                ncdrscr[strand][i] += ncdrscr[strand][i - 1];
            }

            score = intrg.Score(s,ii);
            if (score == BadScore) {
                score = 0;
            }
            ingscr[strand][i] = score;
            if (i > 0) {
                ingscr[strand][i] += ingscr[strand][i - 1];
            }
        }

        for (int frame = 0;  frame < 3;  ++frame) {
            for (int i = 0;  i < len;  ++i)
            {
                int codonshift,ii;
                if (strand == Plus) {
                    // left end of codon is shifted by frame bases to left
                    codonshift = (frame + i)%3;
                    ii = i;
                } else {
                    // right end of codon is shifted by frame bases to right
                    codonshift = (frame - i)%3;
                    if (codonshift < 0) {
                        codonshift += 3;
                    }
                    ii = len - 1-i;
                }

                double score = cdr.Score(s,ii,codonshift);
                if (score == BadScore) {
                    score = 0;
                }

                cdrscr[strand][frame][i] = score;
                if (i > 0) {
                    cdrscr[strand][frame][i] += cdrscr[strand][frame][i - 1];

                    IVec& lstp = laststop[strand][frame];
                    lstp[i] = max(lstp[i - 1],lstp[i]);
                }
            }
        }
    }

    for (int strand = 0;  strand < 2;  ++strand) {
        for (int frame = 0;  frame < 3;  ++frame) {
            for (int i = 0;  i < len;  ++i) {
                cdrscr[strand][frame][i] -= ncdrscr[strand][i];
            }
        }
        for (int i = 0;  i < len;  ++i) {
            ingscr[strand][i] -= ncdrscr[strand][i];

            int left, right;
            if (dscr[strand][i] != BadScore) {
                Terminal& t = donor;
                left = i + 1-(strand == Plus ? t.Left():t.Right());
                right = i + (strand == Plus ? t.Right():t.Left());
                dscr[strand][i] -= NonCodingScore(left,right,strand);
            }
            if (ascr[strand][i] != BadScore) {
                Terminal& t = acceptor;
                left = i + 1-(strand == Plus ? t.Left():t.Right());
                right = i + (strand == Plus ? t.Right():t.Left());
                ascr[strand][i] -= NonCodingScore(left,right,strand);
            }
            if (sttscr[strand][i] != BadScore) {
                Terminal& t = start;
                left = i + 1-(strand == Plus ? t.Left():t.Right());
                right = i + (strand == Plus ? t.Right():t.Left());
                sttscr[strand][i] -= NonCodingScore(left,right,strand);
            }
            if (stpscr[strand][i] != BadScore) {
                Terminal& t = stop;
                left = i + 1-(strand == Plus ? t.Left():t.Right());
                right = i + (strand == Plus ? t.Right():t.Left());
                stpscr[strand][i] -= NonCodingScore(left,right,strand);
            }
        }
    }
}


END_NCBI_SCOPE


/*
 * ===========================================================================
 * $Log: models.cpp,v $
 * Revision 1000.2  2004/06/01 18:08:35  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.3
 *
 * Revision 1.3  2004/05/21 21:41:03  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.2  2003/11/06 00:52:38  ucko
 * Don't try to take the log of an integer constant.
 *
 * Revision 1.1  2003/10/24 15:07:25  dicuccio
 * Initial revision
 *
 * ===========================================================================
 */