test_helper.cpp

ncbi源码

/////////////////////////////////////////////////////////////////////////////
/************************************************************************
    1.2.2. Test function
    Input: scope + seq_id
        other parameters - only to verify results.
    Operations:
        investigate SeqMap, verify lengths
        investigate SeqVector, verify seq.data (as strings)
        investigate SeqMap again - some references may be resolved now
        enumerate descriptions
        enumerate features for the whole sequence
        enumerate features for an interval
        enumerate alignments for the whole sequence
        enumerate alignments for an interval
************************************************************************/
#define CHECK_WRAP() \
    {{ \
        bool got_exception = false; \
        try {
#define CHECK_END2(MSG, have_errors) \
        } catch (exception) { \
            got_exception = true; \
            if ( !have_errors ) { \
                LOG_POST("Can not " MSG); throw; \
            } \
        } \
        if ( have_errors && !got_exception ) { \
            THROW1_TRACE(runtime_error, \
                         "Managed to " MSG " of erroneous sequence"); \
        } \
    }}
#define CHECK_END(MSG)         CHECK_END2(MSG, have_errors)
#define CHECK_END_ALWAYS(MSG)  CHECK_END2(MSG, false)


bool CTestHelper::sm_DumpFeatures = false;
#if defined _MT
bool CTestHelper::sm_TestRemoveEntry = false;
#else 
bool CTestHelper::sm_TestRemoveEntry = true;
#endif

void CTestHelper::ProcessBioseq(CScope& scope, CSeq_id& id,
                                TSeqPos seq_len,
                                string seq_str, string seq_str_compl,
                                int seq_desc_cnt,
                                int seq_feat_ra_cnt, // resolve-all method
                                int seq_feat_cnt, int seq_featrg_cnt,
                                int seq_align_cnt, int seq_alignrg_cnt,
                                size_t feat_annots_cnt,
                                size_t featrg_annots_cnt,
                                size_t align_annots_cnt,
                                size_t alignrg_annots_cnt,
                                bool tse_feat_test,
                                bool have_errors)
{
    CBioseq_Handle handle = scope.GetBioseqHandle(id);
    if ( !handle ) {
        LOG_POST("No seq-id found");
        return;
    }
    handle.GetTopLevelSeqEntry();
#if 0 // build order issues
    CHECK_WRAP();
    sequence::GetTitle(handle);
    CHECK_END("get sequence title");
#endif

    CBioseq_Handle::TBioseqCore seq_core = handle.GetBioseqCore();
    CHECK_WRAP();
    CConstRef<CSeqMap> seq_map(&handle.GetSeqMap());
    TSeqPos len = 0;
    _TRACE("ProcessBioseq("<<id.AsFastaString()<<") seq_len="<<seq_len<<"):");
    // Iterate seq-map except the last element
    len = 0;
    CSeqMap::const_iterator seg = seq_map->begin(&scope);
    vector<CSeqMap::const_iterator> itrs;
    for ( ; seg != seq_map->end(&scope); ++seg ) {
        _ASSERT(seg);
        itrs.push_back(seg);
        switch (seg.GetType()) {
        case CSeqMap::eSeqData:
            _TRACE('@'<<len<<": seqData("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqRef:
            _TRACE('@'<<len<<": seqRef("<<seg.GetLength()<<
                   ", id="<<seg.GetRefSeqid().AsString()<<
                   ", pos="<<seg.GetRefPosition()<<
                   ", minus="<<seg.GetRefMinusStrand()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqGap:
            _TRACE('@'<<len<<": seqGap("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqEnd:
            _ASSERT("Unexpected END segment" && 0);
            break;
        default:
            _ASSERT("Unexpected segment type" && 0);
            break;
        }
    }
    _ASSERT(!seg);
    for ( int i = itrs.size()-1; i >= 0; --i ) {
        --seg;
        _ASSERT(seg);
        _ASSERT(seg == itrs[i]);
    }
    --seg;
    _ASSERT(!seg);

    _TRACE("ProcessBioseq("<<id.AsFastaString()<<") len="<<len<<")");
    _ASSERT(len == seq_len);
    CHECK_END("get sequence map");

    CHECK_WRAP();
    CConstRef<CSeqMap> seq_map(&handle.GetSeqMap());
    TSeqPos len = 0;
    _TRACE("ProcessBioseq("<<id.AsFastaString()<<") - depth-restricted");
    // Iterate seq-map except the last element
    len = 0;
    CSeqMap::const_iterator seg(seq_map, &scope, SSeqMapSelector()
                                .SetResolveCount(2));
    vector<CSeqMap::const_iterator> itrs;
    for ( ; seg != seq_map->end(&scope); ++seg ) {
        _ASSERT(seg);
        itrs.push_back(seg);
        switch (seg.GetType()) {
        case CSeqMap::eSeqData:
            _TRACE('@'<<len<<": seqData("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqRef:
            _TRACE('@'<<len<<": seqRef("<<seg.GetLength()<<
                   ", id="<<seg.GetRefSeqid().AsString()<<
                   ", pos="<<seg.GetRefPosition()<<
                   ", minus="<<seg.GetRefMinusStrand()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqGap:
            _TRACE('@'<<len<<": seqGap("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqEnd:
            _ASSERT("Unexpected END segment" && 0);
            break;
        default:
            _ASSERT("Unexpected segment type" && 0);
            break;
        }
    }
    _ASSERT(!seg);
    for ( int i = itrs.size()-1; i >= 0; --i ) {
        --seg;
        _ASSERT(seg);
        _ASSERT(seg == itrs[i]);
    }
    --seg;
    _ASSERT(!seg);
    CHECK_END("get restricted sequence map");

    CHECK_WRAP();
    CConstRef<CSeqMap> seq_map(&handle.GetSeqMap());
    TSeqPos len = 0;
    _TRACE("ProcessBioseq("<<id.AsFastaString()<<") - TSE-restricted");
    // Iterate seq-map except the last element
    len = 0;
    CSeqMap::const_iterator seg(seq_map, &scope, SSeqMapSelector()
                                .SetLimitTSE(&handle.GetTopLevelSeqEntry())
                                .SetResolveCount(kInvalidSeqPos)
                                .SetFlags(CSeqMap::fFindAny));
    vector<CSeqMap::const_iterator> itrs;
    for ( ; seg != seq_map->end(&scope); ++seg ) {
        _ASSERT(seg);
        itrs.push_back(seg);
        switch (seg.GetType()) {
        case CSeqMap::eSeqData:
            _TRACE('@'<<len<<": seqData("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqRef:
            _TRACE('@'<<len<<": seqRef("<<seg.GetLength()<<
                   ", id="<<seg.GetRefSeqid().AsString()<<
                   ", pos="<<seg.GetRefPosition()<<
                   ", minus="<<seg.GetRefMinusStrand()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqGap:
            _TRACE('@'<<len<<": seqGap("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqEnd:
            _ASSERT("Unexpected END segment" && 0);
            break;
        default:
            _ASSERT("Unexpected segment type" && 0);
            break;
        }
    }
    _ASSERT(!seg);
    for ( int i = itrs.size()-1; i >= 0; --i ) {
        --seg;
        _ASSERT(seg);
        _ASSERT(seg == itrs[i]);
    }
    --seg;
    _ASSERT(!seg);
    CHECK_END("get restricted sequence map");

    CHECK_WRAP();
    CConstRef<CSeqMap> seq_map(&handle.GetSeqMap());
    TSeqPos len = 0;
    _TRACE("ProcessBioseq("<<id.AsFastaString()<<
           ") seq_len="<<seq_len<<") resolved:");
    // Iterate seq-map except the last element
    len = 0;
    CSeqMap::const_iterator seg = seq_map->begin_resolved(&scope);
    vector<CSeqMap::const_iterator> itrs;
    for ( ; seg != seq_map->end_resolved(&scope); ++seg ) {
        _ASSERT(seg);
        itrs.push_back(seg);
        switch (seg.GetType()) {
        case CSeqMap::eSeqData:
            _TRACE('@'<<len<<": seqData("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqRef:
            _TRACE('@'<<len<<": seqRef("<<seg.GetLength()<<
                   ", id="<<seg.GetRefSeqid().AsString()<<
                   ", pos="<<seg.GetRefPosition()<<
                   ", minus="<<seg.GetRefMinusStrand()<<")");
            _ASSERT("Unexpected REF segment" && 0);
            //len += seg.GetLength();
            break;
        case CSeqMap::eSeqGap:
            _TRACE('@'<<len<<": seqGap("<<seg.GetLength()<<")");
            len += seg.GetLength();
            break;
        case CSeqMap::eSeqEnd:
            _ASSERT("Unexpected END segment" && 0);
            break;
        default:
            _ASSERT("Unexpected segment type" && 0);
            break;
        }
    }
    _ASSERT(!seg);
    for ( int i = itrs.size()-1; i >= 0; --i ) {
        --seg;
        _ASSERT(seg);
        _ASSERT(seg == itrs[i]);
    }
    --seg;
    _ASSERT(!seg);

    _TRACE("ProcessBioseq("<<id.AsFastaString()<<") len="<<len<<")");
    _ASSERT(len == seq_len);
    CHECK_END("get resolved sequence map");

    CSeqVector seq_vect;

    CHECK_WRAP();
    seq_vect = handle.GetSeqVector(CBioseq_Handle::eCoding_NotSet);
    string sout = "";
    {{
        CSeqVector_CI vit(seq_vect, 0);
        for ( ; vit.GetPos() < seq_vect.size(); ++vit) {
            sout += *vit;
        }
        _ASSERT(NStr::PrintableString(sout) == seq_str);
        seq_vect.SetCoding(CBioseq_Handle::eCoding_Ncbi);
        int pos = seq_vect.size() - 1;
        vit.SetPos(pos);
        for ( ; pos > 0; --vit, --pos) {
            _ASSERT(sout[vit.GetPos()] == *vit);
        }
        _ASSERT(sout[vit.GetPos()] == *vit);
        vit.GetSeqData(0, seq_vect.size(), sout);
        _ASSERT(NStr::PrintableString(sout) == seq_str);
        sout = "";
        seq_vect.SetCoding(CBioseq_Handle::eCoding_NotSet);
    }}
    for (TSeqPos i = 0; i < seq_vect.size(); i++) {
        sout += seq_vect[i];
    }
    _ASSERT(NStr::PrintableString(sout) == seq_str);
    CHECK_END("get seq vector");

    CHECK_WRAP();
    CSeq_interval interval;
    interval.SetId(id);
    interval.SetFrom(seq_str.size());
    interval.SetTo(seq_str.size()-1);
    CSeq_loc loc;
    loc.SetInt(interval);
    seq_vect = handle.GetSequenceView(loc, CBioseq_Handle::eViewConstructed);
    string sout = "";
    for (TSeqPos i = 0; i < seq_vect.size(); i++) {
        sout += seq_vect[i];
    }
    _ASSERT(sout.empty());
    CHECK_END_ALWAYS("get seq view");

    if (seq_core->GetInst().IsSetStrand() &&
        seq_core->GetInst().GetStrand() == CSeq_inst::eStrand_ds) {
        CHECK_WRAP();
        seq_vect = handle.GetSeqVector(CBioseq_Handle::eCoding_NotSet,
                                       eNa_strand_minus);
        string sout_rev = "";
        for (TSeqPos i = seq_vect.size(); i> 0; i--) {
            sout_rev += seq_vect[i-1];
        }
        _ASSERT(NStr::PrintableString(sout_rev) == seq_str_compl);
        CHECK_END("get reverse seq vector");
    }
    else {
        _ASSERT(seq_str_compl.empty());
    }

    CConstRef<CBioseq> bioseq(&(handle.GetBioseq()));

    int count = 0;
    CHECK_WRAP();
    // Test CSeq_descr iterator
    for (CSeq_descr_CI desc_it(handle); desc_it;  ++desc_it) {
        count++;
        //### _ASSERT(desc_it->
    }
    _ASSERT(count == seq_desc_cnt);
    CHECK_END_ALWAYS("count CSeq_descr");

    CHECK_WRAP();
    // Test CSeq_feat iterator, resolve-all method
    CSeq_loc loc;
    loc.SetWhole(id);
    count = 0;
    if ( !tse_feat_test ) {
        if ( sm_DumpFeatures ) {
            NcbiCout << "-------------------- CFeat_CI resolve-all over "<<
                id.AsFastaString()<<" --------------------" << NcbiEndl;
        }
        for ( CFeat_CI feat_it(scope, loc, SAnnotSelector().SetResolveAll());
              feat_it;  ++feat_it) {
            count++;
            if ( sm_DumpFeatures ) {
                auto_ptr<CObjectOStream>
                    out(CObjectOStream::Open(eSerial_AsnText, NcbiCout));
                *out << feat_it->GetMappedFeature();
            }
            //### _ASSERT(feat_it->
        }
        _ASSERT(count == seq_feat_ra_cnt);
    }
    CHECK_END("get annot set");

    CHECK_WRAP();
    // Test CSeq_feat iterator
    CSeq_loc loc;
    loc.SetWhole(id);
    count = 0;
    set<const CSeq_annot*> annot_set;
    if ( sm_DumpFeatures ) {
        NcbiCout << "-------------------- CFeat_CI over "<<
            id.AsFastaString()<<" --------------------" << NcbiEndl;
    }
    if ( !tse_feat_test ) {
        for ( CFeat_CI feat_it(scope, loc); feat_it;  ++feat_it) {
            count++;
            annot_set.insert(&feat_it.GetSeq_annot());
            if ( sm_DumpFeatures ) {
                auto_ptr<CObjectOStream>
                    out(CObjectOStream::Open(eSerial_AsnText, NcbiCout));
                *out << feat_it->GetMappedFeature();
            }
            //### _ASSERT(feat_it->
        }
        if ( sm_DumpFeatures ) {
            NcbiCout << "-------------------- "
                "product CFeat_CI --------------------" << NcbiEndl;
        }
        // Get products
        for ( CFeat_CI feat_it(scope, loc, SAnnotSelector().SetByProduct());
              feat_it;  ++feat_it) {
            count++;
            annot_set.insert(&feat_it.GetSeq_annot());
            if ( sm_DumpFeatures ) {
                auto_ptr<CObjectOStream>
                    out(CObjectOStream::Open(eSerial_AsnText, NcbiCout));
                *out << feat_it->GetMappedFeature();
            }
            //### _ASSERT(feat_it->
        }
    }
    else {
        if ( sm_DumpFeatures ) {
            NcbiCout << "-------------------- "
                "LimitTSE --------------------" << NcbiEndl;
        }
        for ( CFeat_CI feat_it(handle, 0, 0,
                               SAnnotSelector()
                               .SetLimitTSE(&handle.GetTopLevelSeqEntry()));
              feat_it;  ++feat_it) {
            count++;
            annot_set.insert(&feat_it.GetSeq_annot());
            if ( sm_DumpFeatures ) {