seq_align.cpp

ncbi源码

    if ( !GetSegs().IsStd() ) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateDensegFromStdseg(): "
                   "Input Seq-align should have segs of type StdSeg!");
    }

    CRef<CSeq_align> sa(new CSeq_align);
    sa->SetType(eType_not_set);
    CDense_seg& ds = sa->SetSegs().SetDenseg();

    typedef CDense_seg::TDim    TNumrow;
    typedef CDense_seg::TNumseg TNumseg;

    vector<TSeqPos>       row_lens;
    CDense_seg::TLens&    seg_lens = ds.SetLens();
    CDense_seg::TStarts&  starts   = ds.SetStarts();
    CDense_seg::TStrands& strands  = ds.SetStrands();
    CDense_seg::TWidths&  widths   = ds.SetWidths();
    vector<bool>          widths_determined;

    TSeqPos row_len;
    TSeqPos from, to;
    ENa_strand strand;


    TNumseg seg = 0;
    TNumrow dim = 0, row = 0;
    ITERATE (C_Segs::TStd, std_i, GetSegs().GetStd()) {

        const CStd_seg& ss = **std_i;

        seg_lens.push_back(0);
        TSeqPos& seg_len = seg_lens.back();
        row_len = 0;
        row_lens.clear();
        widths_determined.push_back(false);

        row = 0;
        ITERATE (CStd_seg::TLoc, i, ss.GetLoc()) {

            const CSeq_id* seq_id = 0;
            // push back initialization values
            if (seg == 0) {
                widths.push_back(0);
                strands.push_back(eNa_strand_unknown);
            }

            if ((*i)->IsInt()) {
                const CSeq_interval& interval = (*i)->GetInt();
                
                // determine start and len
                from = interval.GetFrom();
                to = interval.GetTo();
                starts.push_back(from);
                row_len = to - from + 1;
                row_lens.push_back(row_len);
                
                int width = 0;
                // try to determine/check the seg_len and width
                if (!seg_len) {
                    width = 0;
                    seg_len = row_len;
                } else {
                    if (row_len * 3 == seg_len) {
                        seg_len /= 3;
                        width = 1;
                    } else if (row_len / 3 == seg_len) {
                        width = 3;
                    } else if (row_len != seg_len) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    }
                }
                if (width > 0) {
                    widths_determined[seg] = true;
                    if (widths[row] > 0  &&  widths[row] != width) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    } else {
                        widths[row] = width;
                    }
                }

                // get the id
                seq_id = &(*i)->GetInt().GetId();

                // determine/check the strand
                if (interval.CanGetStrand()) {
                    strand = interval.GetStrand();
                    if (seg == 0  ||  strands[row] == eNa_strand_unknown) {
                        strands[row] = strand;
                    } else {
                        if (strands[row] != strand) {
                            NCBI_THROW(CSeqalignException,
                                       eInvalidInputAlignment,
                                       "CreateDensegFromStdseg(): "
                                       "Inconsistent strands!");
                        }
                    }
                } else {
                    strand = eNa_strand_unknown;
                }

                    
            } else if ((*i)->IsEmpty()) {
                starts.push_back(-1);
                if (seg == 0) {
                    strands[row] = eNa_strand_unknown;
                }
                seq_id = &(*i)->GetEmpty();
                row_lens.push_back(0);

            }

            // determine/check the id
            if (seg == 0) {
                CRef<CSeq_id> id(new CSeq_id);
                SerialAssign(*id, *seq_id);
                ds.SetIds().push_back(id);
            } else {
                CSeq_id& id = *ds.SetIds()[row];
                if (!SerialEquals(id, *seq_id)) {
                    if (SeqIdChooser) {
                        SeqIdChooser->ChooseSeqId(id, *seq_id);
                    } else {
                        string errstr =
                            string("CreateDensegFromStdseg(): Seq-ids: ") +
                            id.AsFastaString() + " and " +
                            seq_id->AsFastaString() + " are not identical!" +
                            " Without the OM it cannot be determined if they belong to"
                            " the same sequence."
                            " Define and pass ChooseSeqId to resolve seq-ids.";
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment, errstr);
                    }
                }
            }

            // next row
            row++;
            if (seg == 0) {
                dim++;
            }
        }
        if (dim != ss.GetDim()  ||  row != dim) {
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                       "CreateDensegFromStdseg(): "
                       "Inconsistent dimentions!");
        }

        if (widths_determined[seg]) {
            // go back and determine/check widths
            for (row = 0; row < dim; row++) {
                if ((row_len = row_lens[row]) > 0) {
                    int width = 0;
                    if (row_len == seg_len * 3) {
                        width = 3;
                    } else if (row_len == seg_len) {
                        width = 1;
                    }
                    if (widths[row] > 0  &&  widths[row] != width) {
                        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                                   "CreateDensegFromStdseg(): "
                                   "Std-seg segment lengths not accurate!");
                    } else {
                        widths[row] = width;
                    }
                }
            }
        }

        // next seg
        seg++;
    }

    ds.SetDim(dim);
    ds.SetNumseg(seg);

    // go back and finish lens determination
    bool widths_failure = false;
    bool widths_success = false;
    for (seg = 0; seg < ds.GetNumseg(); seg++) {
        if (!widths_determined[seg]) {
            for(row = 0; row < dim; row++) {
                if (starts[seg * dim + row] >= 0) {
                    int width = widths[row];
                    if (width == 3) {
                        seg_lens[seg] /= 3;
                    } else if (width == 0) {
                        widths_failure = true;
                    }
                    break;
                }
            }
        } else {
            widths_success = true;
        }
    }

    if (widths_failure) {
        if (widths_success) {
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                       "CreateDensegFromStdseg(): "
                       "Some widths cannot be determined!");
        } else {
            ds.ResetWidths();
        }
    }
    
    // finish the strands
    for (seg = 1; seg < ds.GetNumseg(); seg++) {
        for (row = 0; row < dim; row++) {
            strands.push_back(strands[row]);
        }
    }

    return sa;
}


///----------------------------------------------------------------------------
/// PRE : the Seq-align is a Dense-seg of aligned nucleotide sequences
/// POST: Seq_align of type Dense-seg is created with each of the m_Widths 
///       equal to 3 and m_Lengths devided by 3.
CRef<CSeq_align> 
CSeq_align::CreateTranslatedDensegFromNADenseg() const
{
    if ( !GetSegs().IsDenseg() ) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateTranslatedDensegFromNADenseg(): "
                   "Input Seq-align should have segs of type Dense-seg!");
    }
    
    CRef<CSeq_align> sa(new CSeq_align);
    sa->SetType(eType_not_set);

    if (GetSegs().GetDenseg().IsSetWidths()) {
        NCBI_THROW(CSeqalignException, eInvalidInputAlignment,
                   "CSeq_align::CreateTranslatedDensegFromNADenseg(): "
                   "Widths already exist for the original alignment");
    }

    // copy from the original
    sa->Assign(*this);

    CDense_seg& ds = sa->SetSegs().SetDenseg();

    // fix the lengths
    const CDense_seg::TLens& orig_lens = GetSegs().GetDenseg().GetLens();
    CDense_seg::TLens&       lens      = ds.SetLens();

    for (CDense_seg::TNumseg numseg = 0; numseg < ds.GetNumseg(); numseg++) {
        if (orig_lens[numseg] % 3) {
            string errstr =
                string("CSeq_align::CreateTranslatedDensegFromNADenseg(): ") +
                "Length of segment " + NStr::IntToString(numseg) +
                " is not divisible by 3.";
            NCBI_THROW(CSeqalignException, eInvalidInputAlignment, errstr);
        } else {
            lens[numseg] = orig_lens[numseg] / 3;
        }
    }

    // add the widths
    ds.SetWidths().resize(ds.GetDim(), 3);

#if _DEBUG
    ds.Validate(true);
#endif

    return sa;
}


END_objects_SCOPE // namespace ncbi::objects::

END_NCBI_SCOPE


/*
* ===========================================================================
*
* $Log: Seq_align.cpp,v $
* Revision 1000.2  2004/06/01 19:33:39  gouriano
* PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.16
*
* Revision 1.16  2004/05/19 17:25:43  gorelenk
* Added include of PCH - ncbi_pch.hpp
*
* Revision 1.15  2004/05/05 19:16:25  johnson
* Added SwapRows method for 'disc' seq-align / seq-align-set
*
* Revision 1.14  2004/04/27 19:17:13  johnson
* Added GetNamedScore helper function
*
* Revision 1.13  2004/04/19 17:27:22  grichenk
* Added GetSeq_id(TDim row)
*
* Revision 1.12  2004/04/12 16:51:02  vasilche
* Fixed uninitialized variables.
*
* Revision 1.11  2004/03/09 17:27:42  kuznets
* Bug fix(compilation): CSeq_align::CreateDensegFromStdseg MSVC
*
* Revision 1.10  2004/03/09 17:14:33  todorov
* changed the C-style callback to a functor
*
* Revision 1.9  2004/02/23 15:31:24  todorov
* +TChooseSeqIdCallback to abstract resolving seq-ids in CreateDensegFromStdseg()
*
* Revision 1.8  2004/02/13 17:10:24  todorov
* - inconsistent ids exception in CreateDensegFromStdseg
*
* Revision 1.7  2004/01/09 16:32:06  todorov
* +SetType(eType_not_set)
*
* Revision 1.6  2003/12/16 22:54:31  todorov
* +CreateTranslatedDensegFromNADenseg
*
* Revision 1.5  2003/12/11 22:34:05  todorov
* Implementation of GetSeq{Start,Stop}
*
* Revision 1.4  2003/09/16 15:31:14  todorov
* Added validation methods. Added seq range methods
*
* Revision 1.3  2003/08/26 20:28:38  johnson
* added 'SwapRows' method
*
* Revision 1.2  2003/08/19 21:11:13  todorov
* +CreateDensegFromStdseg
*
* Revision 1.1  2003/08/13 18:12:03  johnson
* added 'Reverse' method
*
*
* ===========================================================================
*/
/* Original file checksum: lines: 64, chars: 1885, CRC32: 4e5d1825 */