alnmix.cpp

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

                }
            }
        }
    }
    x_CreateRowsVector();
    x_CreateSegmentsVector();
    x_CreateDenseg();
}

CAlnMix::TSecondRowFits
CAlnMix::x_SecondRowFits(CAlnMixMatch * match) const
{
    CAlnMixSeq::TStarts&          starts1 = match->m_AlnSeq1->m_Starts;
    CAlnMixSeq::TStarts&          starts2 = match->m_AlnSeq2->m_Starts;
    CAlnMixSeq                  * seq1    = match->m_AlnSeq1,
                                * seq2    = match->m_AlnSeq2;
    TSeqPos&                      start1  = match->m_Start1;
    TSeqPos&                      start2  = match->m_Start2;
    TSeqPos&                      len     = match->m_Len;
    const int&                    width1  = seq1->m_Width;
    const int&                    width2  = seq2->m_Width;
    CAlnMixSeq::TStarts::iterator start_i;
    TSignedSeqPos                 delta, delta1, delta2;

    // subject sequences go on separate rows if requested
    if (m_MergeFlags & fQuerySeqMergeOnly) {
        if (seq2->m_DSIndex) {
            if (seq2->m_DSIndex == match->m_DSIndex) {
                return eSecondRowFitsOk;
            } else {
                return eForceSeparateRow;
            }
        } else {
            seq2->m_DSIndex = match->m_DSIndex;
            return eSecondRowFitsOk;
        }
    }

    if ( !starts2.empty() ) {

        // check strand
        if (seq2->m_PositiveStrand !=
            (seq1->m_PositiveStrand ?
             !match->m_StrandsDiffer :
             match->m_StrandsDiffer)) {
            return eInconsistentStrand;
        }

        // check frame
        if (seq2->m_Width == 3  &&  seq2->m_Frame != start2 % 3) {
            return eInconsistentFrame;
        }

        start_i = starts2.lower_bound(start2);

        // check below
        if (start_i != starts2.begin()) {
            start_i--;
            
            // check for inconsistency on the first row
            if ( !m_IndependentDSs ) {
                CAlnMixSegment::TStartIterators::iterator start_it_i =
                    start_i->second->m_StartIts.find(seq1);
                if (start_it_i != start_i->second->m_StartIts.end()) {
                    if (match->m_StrandsDiffer) {
                        delta = start1 + len * width1 - start_it_i->second->first;
                        if (delta > 0) {
                            // target above
                            // x----- x-)-)
                            // (----x (---x
                            // target below
                            if (m_MergeFlags & fTruncateOverlaps) {
                                delta /= width1;
                                if ((len -= delta) > 0) {
                                    start2 += delta * width2;
                                } else {
                                    return eIgnoreMatch;
                                }
                            } else {
                                return eFirstRowOverlapAbove;
                            }
                        }
                    } else {
                        delta = start_it_i->second->first
                            + start_i->second->m_Len * width1
                            - start1;
                        if (delta > 0) {
                            // below target
                            // x---- x-)--)
                            // x---) x----)
                            // below target
                            if (m_MergeFlags & fTruncateOverlaps) {
                                delta /= width1;
                                if ((len -= delta) > 0) {
                                    start1 += delta * width1;
                                    start2 += delta * width2;
                                } else {
                                    return eIgnoreMatch;
                                }
                            } else {
                                return eFirstRowOverlapBelow;
                            }
                        }
                    }
                }
            }

            // check for overlap with the segment below on second row
            if ((delta = start_i->first + start_i->second->m_Len * width2
                 - start2) > 0) {
                //       target
                // ----- ------
                // x---- x-)--)
                // below target
                if (m_MergeFlags & fTruncateOverlaps) {
                    delta /= width2;
                    if ((len -= delta) > 0) {
                        start2 += delta * width2;
                        if ( !match->m_StrandsDiffer ) {
                            start1 += delta * width1;
                        }
                    } else {
                        return eIgnoreMatch;
                    }
                } else {
                    return eSecondRowOverlap;
                }
            }
            if (start_i != starts2.end()) {
                start_i++;
            }
        }

        // check the overlapping area for consistency
        while (start_i != starts2.end()  &&  
               start_i->first < start2 + len * width2) {
            if ( !m_IndependentDSs ) {
                CAlnMixSegment::TStartIterators::iterator start_it_i =
                    start_i->second->m_StartIts.find(seq1);
                if (start_it_i != start_i->second->m_StartIts.end()) {
                    if (match->m_StrandsDiffer) {
                        // x---..- x---..--)
                        // (---..- (--x..--x
                        delta1 = (start1 - start_it_i->second->first) / width1 +
                            len - start_i->second->m_Len;
                    } else {
                        // x--..- x---...-)
                        // x--..- x---...-)
                        delta1 = (start_it_i->second->first - start1) / width1;
                    }
                    delta2 = (start_i->first - start2) / width2;
                    if (delta1 != delta2) {
                        if (m_MergeFlags & fTruncateOverlaps) {
                            delta = (delta1 < delta2 ? delta1 : delta2);
                            if (match->m_StrandsDiffer) {
                                start1 += (len - delta) * width1;
                            }
                            len = delta;
                        } else {
                            return eInconsistentOverlap;
                        }
                    }
                }
            }
            start_i++;
        }

        // check above for consistency
        if (start_i != starts2.end()) {
            if ( !m_IndependentDSs ) {
                CAlnMixSegment::TStartIterators::iterator start_it_i =
                    start_i->second->m_StartIts.find(seq1);
                if (start_it_i != start_i->second->m_StartIts.end()) {
                    if (match->m_StrandsDiffer) {
                        delta = start_it_i->second->first + 
                            start_i->second->m_Len * width1 - start1;
                        if (delta > 0) {
                            // below target
                            // x---- x-)--)
                            // (---x (----x
                            // above target
                            if (m_MergeFlags & fTruncateOverlaps) {
                                if ((len -= delta / width1) > 0) {
                                    start1 += delta;
                                } else {
                                    return eIgnoreMatch;
                                }
                            } else {
                                return eFirstRowOverlapBelow;
                            }
                        }
                    } else {
                        delta = start1 + len * width1 - start_it_i->second->first;
                        if (delta > 0) {
                            // target above
                            // x--x-) ----)
                            // x----) x---)
                            // target above
                            if (m_MergeFlags & fTruncateOverlaps) {
                                if ((len -= delta / width1) <= 0) {
                                    return eIgnoreMatch;
                                }
                            } else {
                                return eFirstRowOverlapAbove;
                            }
                        }
                    }
                }
            }
        }

        // check for inconsistent matches
        if ((start_i = starts1.find(start1)) == starts1.end() ||
            start_i->first != start1) {
            // commenting out for now, since moved the function call ahead            
//             NCBI_THROW(CAlnException, eMergeFailure,
//                        "CAlnMix::x_SecondRowFits(): "
//                        "Internal error: starts1 do not match");
        } else {
            CAlnMixSegment::TStartIterators::iterator it;
            TSeqPos tmp_start =
                match->m_StrandsDiffer ? start2 + len * width2 : start2;
            while (start_i != starts1.end()  &&
                   start_i->first < start1 + len * width1) {

                CAlnMixSegment::TStartIterators& its = 
                    start_i->second->m_StartIts;

                if (match->m_StrandsDiffer) {
                    tmp_start -= start_i->second->m_Len * width2;
                }

                if ((it = its.find(seq2)) != its.end()) {
                    if (it->second->first != tmp_start) {
                        // found an inconsistent prev match
                        return eSecondRowInconsistency;
                    }
                }

                if ( !match->m_StrandsDiffer ) {
                    tmp_start += start_i->second->m_Len * width2;
                }

                start_i++;
            }
        }
    }
    return eSecondRowFitsOk;
}

void CAlnMix::x_CreateRowsVector()
{
    m_Rows.clear();

    int count = 0;
    ITERATE (TSeqs, i, m_Seqs) {
        CAlnMixSeq * seq = *i;
        m_Rows.push_back(seq);
        seq->m_RowIndex = count++;
        while ( (seq = seq->m_ExtraRow) != NULL ) {
            seq->m_RowIndex = count++;
            m_Rows.push_back(seq);
        }
    }
}


void CAlnMix::x_CreateSegmentsVector()
{
    TSegmentsContainer gapped_segs;

    // init the start iterators for each row
    NON_CONST_ITERATE (TSeqs, row_i, m_Rows) {
        CAlnMixSeq * row = *row_i;
        if ( !row->m_Starts.empty() ) {
            if (row->m_PositiveStrand) {
                row->m_StartIt = row->m_Starts.begin();
            } else {
                row->m_StartIt = row->m_Starts.end();
                row->m_StartIt--;
            }
        } else {
            row->m_StartIt = row->m_Starts.end();
#if _DEBUG
            string errstr =
                string("CAlnMix::x_CreateSegmentsVector():") +
                " Internal error: no starts for row " +
                NStr::IntToString(row->m_RowIndex) +
                " (seq " +
                NStr::IntToString(row->m_SeqIndex) + ").";
            NCBI_THROW(CAlnException, eMergeFailure, errstr);
#endif
        }
    }

    // init the start iterators for each extra row
    NON_CONST_ITERATE (list<CRef<CAlnMixSeq> >, row_i, m_ExtraRows) {
        CAlnMixSeq * row = *row_i;
        if ( !row->m_Starts.empty() ) {
            if (row->m_PositiveStrand) {
                row->m_StartIt = row->m_Starts.begin();
            } else {
                row->m_StartIt = row->m_Starts.end();
                row->m_StartIt--;
            }
        } else {
            row->m_StartIt = row->m_Starts.end();
#if _DEBUG
            string errstr =
                string("CAlnMix::x_CreateSegmentsVector():") +
                " Internal error: no starts for row " +
                NStr::IntToString(row->m_RowIndex) + ".";
            NCBI_THROW(CAlnException, eMergeFailure, errstr);
#endif
        }
    }

#if _DEBUG
    ITERATE (TSeqs, row_i, m_Rows) {
        ITERATE (CAlnMixSegment::TStarts, st_i, (*row_i)->m_Starts) {
            ITERATE(CAlnMixSegment::TStartIterators,
                    st_it_i, (*st_i).second->m_StartIts) {
                // both should point to the same seg
                if ((*st_it_i).second->second != (*st_i).second) {
                    string errstr =
                        string("CAlnMix::x_CreateSegmentsVector():")
                        + " Internal error: Segments messed up."
                        + " row=" + NStr::IntToString((*row_i)->m_RowIndex)
                        + " seq=" + NStr::IntToString((*row_i)->m_SeqIndex)
                        + " strand=" +
                        ((*row_i)->m_PositiveStrand ? "plus" : "minus");
                    NCBI_THROW(CAlnException, eMergeFailure, errstr);
                }
            }
        }