dense_seg.cpp

ncbi源码

//-----------------------------------------------------------------------------
// PRE : numbers of the rows to swap
// POST: alignment rearranged with row1 where row2 used to be & vice versa
void CDense_seg::SwapRows(TDim row1, TDim row2)
{
    if (row1 >= GetDim()  ||  row1 < 0  ||
        row2 >= GetDim()  ||  row2 < 0) {
        NCBI_THROW(CSeqalignException, eOutOfRange,
                   "Row numbers supplied to CDense_seg::SwapRows must be "
                   "in the range [0, dim)");
    }

    //swap ids
    swap(SetIds()[row1], SetIds()[row2]);

    int idxStop = GetNumseg()*GetDim();
    
    //swap starts
    for(int i = 0; i < idxStop; i += GetDim()) {
        swap(SetStarts()[i+row1], SetStarts()[i+row2]);
    }

    //swap strands
    for(int i = 0; i < idxStop; i += GetDim()) {
        swap(SetStrands()[i+row1], SetStrands()[i+row2]);
    }
}


void CDense_seg::RemapToLoc(TDim row, const CSeq_loc& loc,
                            bool ignore_strand)
{
    if (loc.IsWhole()) {
        return;
    }

    TSeqPos row_stop  = GetSeqStop(row);

    size_t  ttl_loc_len = 0;
    {{
        CSeq_loc_CI seq_loc_i(loc);
        do {
            ttl_loc_len += seq_loc_i.GetRange().GetLength();
        } while (++seq_loc_i);
    }}

    // check the validity of the seq-loc
    if (ttl_loc_len < row_stop + 1) {
        string errstr("CDense_seg::RemapToLoc():"
                      " Seq-loc is not long enough to"
                      " cover the alignment!"
                      " Maximum row seq pos is ");
        errstr += NStr::IntToString(row_stop);
        errstr += " The total seq-loc len is only ";
        errstr += NStr::IntToString(ttl_loc_len);
        errstr += ", it should be at least ";
        errstr += NStr::IntToString(row_stop+1);
        errstr += " (= max seq pos + 1).";
        NCBI_THROW(CSeqalignException, eOutOfRange, errstr);
    }

    const CDense_seg::TStarts&  starts  = GetStarts();
    const CDense_seg::TStrands& strands = GetStrands();
    const CDense_seg::TLens&    lens    = GetLens();

    const size_t& numrows = CheckNumRows();
    const size_t& numsegs = CheckNumSegs();

    CSeq_loc_CI seq_loc_i(loc);

    TSeqPos start, loc_len, len, len_so_far;
    start = seq_loc_i.GetRange().GetFrom();
    len = loc_len = seq_loc_i.GetRange().GetLength();
    len_so_far = 0;
    
    bool row_plus = !strands.size() || strands[row] != eNa_strand_minus;
    bool loc_plus = seq_loc_i.GetStrand() != eNa_strand_minus;

    // iterate through segments
    size_t  idx = loc_plus ? row : (numsegs - 1) * numrows + row;
    TNumseg seg = loc_plus ? 0 : numsegs - 1;
    while (loc_plus ? seg < GetNumseg() : seg >= 0) {
        if (starts[idx] == -1) {
            // ignore gaps in our sequence
            if (loc_plus) {
                idx += numrows; seg++;
            } else {
                idx -= numrows; seg--;
            }
            continue;
        }

        // iterate the seq-loc if needed
        if ((loc_plus == row_plus ?
             starts[idx] : ttl_loc_len - starts[idx] - lens[seg])
            > len_so_far + loc_len) {

            if (++seq_loc_i) {
                len_so_far += len;
                len   = seq_loc_i.GetRange().GetLength();
                start = seq_loc_i.GetRange().GetFrom();
            } else {
                NCBI_THROW(CSeqalignException, eInvalidInputData,
                           "CDense_seg::RemapToLoc():"
                           " Internal error");
            }

            // assert the strand is the same
            if (loc_plus != (seq_loc_i.GetStrand() != eNa_strand_minus)) {
                NCBI_THROW(CSeqalignException, eInvalidInputData,
                           "CDense_seg::RemapToLoc():"
                           " The strand should be the same accross"
                           " the input seq-loc");
            }
        }

        // offset for the starting position
        if (loc_plus == row_plus) {
            SetStarts()[idx] += start - len_so_far;
        } else {
            SetStarts()[idx] = 
                start - len_so_far + ttl_loc_len - starts[idx] - lens[seg];
        }

        if (lens[seg] > len) {
            TSignedSeqPos len_diff = lens[seg] - len;
            while (1) {
                // move to the next loc part that extends beyond our length
                ++seq_loc_i;
                if (seq_loc_i) {
                    start = seq_loc_i.GetRange().GetFrom();
                } else {
                    NCBI_THROW(CSeqalignException, eOutOfRange,
                               "CDense_seg::RemapToLoc():"
                               " Internal error");
                }

                // split our segment
                SetLens().insert(SetLens().begin() + 
                                 (loc_plus ? seg : seg + 1),
                                 len);
                SetLens()[loc_plus ? seg + 1 : seg] = len_diff;

                // insert new data to account for our split segment
                TStarts temp_starts(numrows, -1);
                for (size_t row_i = 0, tmp_idx = seg * numrows;
                     row_i < numrows;  ++row_i, ++tmp_idx) {
                    TSignedSeqPos& this_start = SetStarts()[tmp_idx];
                    if (this_start != -1) {
                        temp_starts[row_i] = this_start;
                        if (loc_plus == (strands[row_i] != eNa_strand_minus)) {
                            if ((size_t) row == row_i) {
                                temp_starts[row_i] = start;
                            } else {
                                temp_starts[row_i] += len;
                            }
                        } else {
                            this_start += len_diff;
                        }
                    }
                }

                len_so_far += loc_len;
                len = loc_len = seq_loc_i.GetRange().GetLength();

                SetStarts().insert(SetStarts().begin() +
                                   (loc_plus ? seg + 1 : seg) * numrows,
                                   temp_starts.begin(), temp_starts.end());
                
                if (strands.size()) {
                    SetStrands().insert
                        (SetStrands().begin(),
                         strands.begin(), strands.begin() + numrows);
                }

                SetNumseg()++;
                
                if ((len_diff = lens[seg] - len) > 0) {
                    if (loc_plus) {
                        idx += numrows; seg++;
                    } else {
                        idx -= numrows; seg--;
                    }
                } else {
                    break;
                }
            }
        } else {
            len -= lens[seg];
        }

        if (loc_plus) {
            idx += numrows; seg++;
        } else {
            idx -= numrows; seg--;
        }
    } // while iterating through segments
    
    // finally, modify the strands if different
    if ( !ignore_strand ) {
        if (loc_plus != row_plus) {
            if (!strands.size()) {
                // strands do not exist, create them
                SetStrands().resize(GetNumseg() * GetDim(), eNa_strand_plus);
            }
            for (seg = 0, idx = row;
                 seg < GetNumseg(); seg++, idx += numrows) {
                SetStrands()[idx] = loc_plus ? eNa_strand_plus : eNa_strand_minus;
            }
        }
    }

}


END_objects_SCOPE // namespace ncbi::objects::

END_NCBI_SCOPE


/*
* ===========================================================================
*
* $Log: Dense_seg.cpp,v $
* Revision 1000.4  2004/06/01 19:33:35  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.13
*
* Revision 1.13  2004/05/19 17:25:43  gorelenk
* Added include of PCH - ncbi_pch.hpp
*
* Revision 1.12  2004/05/06 18:23:54  todorov
* + optional ignore_strand param to RemapToLoc
*
* Revision 1.11  2004/03/10 13:04:18  dicuccio
* Compilation fix for Win32: don't add char* and char*
*
* Revision 1.10  2004/03/09 21:57:03  todorov
* Fixed the out-of-range exception txt
*
* Revision 1.9  2004/02/12 20:54:15  yazhuk
* Fixed GetStartPos() GetStopPos() handling of empty m_Strands
*
* Revision 1.8  2003/12/19 20:15:21  todorov
* RemapToLoc should do nothing in case of a whole Seq-loc
*
* Revision 1.7  2003/11/20 21:26:33  todorov
* RemapToLoc bug fixes: + seg inc; + loc_len vs len
*
* Revision 1.6  2003/11/04 14:44:46  todorov
* +RemapToLoc
*
* Revision 1.5  2003/09/25 17:50:14  dicuccio
* Changed testing of STL container size to use empty() and avoid warning on MSVC
*
* Revision 1.4  2003/09/16 15:31:14  todorov
* Added validation methods. Added seq range methods
*
* Revision 1.3  2003/08/26 21:10:49  ucko
* #include Seq_id.hpp
*
* Revision 1.2  2003/08/26 20:28:38  johnson
* added 'SwapRows' method
*
* Revision 1.1  2003/08/13 18:12:03  johnson
* added 'Reverse' method
*
*
* ===========================================================================
*/
/* Original file checksum: lines: 64, chars: 1885, CRC32: 4483973b */