utils.cpp

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

        case CSeqFeatData::eSubtype_mRNA:
            parent_type = CSeqFeatData::eSubtype_gene;
            break;

        default:
            // don't link
            out_feats.push_back(*iter);
            continue;
        }

        //
        // find the suitable pool of features we can look at
        //

        bool linked = false;
        CLayoutFeat::TFeatList::reverse_iterator start(iter);
        CLayoutFeat::TFeatList::reverse_iterator end  (feats.begin());
        for ( ;  start != end;  ++start) {
            CLayoutFeat& curr = **start;

            // only link features if they live in the same annotation
            if (curr.GetMappedFeature().GetAnnot().IsNamed()  &&
                feat.GetMappedFeature().GetAnnot().IsNamed()  &&
                curr.GetMappedFeature().GetAnnot().GetName() !=
                feat.GetMappedFeature().GetAnnot().GetName()) {
                continue;
            }

            // make sure we're looking at the right parent type
            CSeqFeatData::ESubtype type =
                curr.GetFeature().GetData().GetSubtype();
            if (type != parent_type) {
                continue;
            }

            // an mRNA can have only one child bound to it
            if (type == CSeqFeatData::eSubtype_mRNA  &&
                curr.GetChildren().size() >= 1) {
                continue;
            }

            // only link features if their locations are swallowable
            sequence::ECompare comp =
                sequence::Compare(curr.GetLocation(), feat.GetLocation());
            if (comp != sequence::eContains  &&  comp != sequence::eSame) {
                continue;
            }

            feat.SetParent(&curr);
            curr.SetChildren().push_back(*iter);
            linked = true;
            break;
        }

        if ( !linked ) {
            out_feats.push_back(*iter);
        }
    }

    out_feats.swap(feats);
}


//
// GetAlignments()
// this retrieves the of alignments in a document by annotation
// as a set of smears.
//
void CSeqUtils::GetAlignmentSmears(const CBioseq_Handle&    handle,
                                   const TSeqRange&         range,
                                   TSeqPos                  window,
                                   CLayoutEngine::TObjects& aligns)
{
    aligns.clear();
    if ( !handle ) {
        return;
    }

    SAnnotSelector sel_annot;
    CAnnot_CI   annot_iter(handle, range.GetFrom(), range.GetTo(), sel_annot);
    
    for (int an_i = 0;  annot_iter;  ++annot_iter, ++an_i)  {
    
        // TODO: remember the annotation's name. That is the smear's label.
        if (annot_iter->IsNamed()) {
            // cout << an_i << " Named: " << annot_iter->GetName() << endl;
        } else {
            // cout << an_i << " Unnamed: " << endl;
        }
        
        const CSeq_annot& seq_annot = annot_iter->GetSeq_annot();
        if ( ! seq_annot.GetData().IsAlign()) {
            continue;
        }
        
        
        CRef<CLayoutObject> ref;

        if (CAlignmentSmear::SeparateStrands(seq_annot)) {
            ref.Reset(new CLayoutAlignSmear(handle,
                                            range.GetFrom(), range.GetTo(), 
                                            window,
                                            CAlignmentSmear::eSmearStrand_Pos,
                                            seq_annot));
            if (ref) {
                aligns.push_back(ref);
            }
            ref.Reset(new CLayoutAlignSmear(handle,
                                            range.GetFrom(), range.GetTo(), 
                                            window,
                                            CAlignmentSmear::eSmearStrand_Neg,
                                            seq_annot));
            if (ref) {
                aligns.push_back(ref);
            }
        } else {
            ref.Reset(new CLayoutAlignSmear(handle,
                                            range.GetFrom(), range.GetTo(), 
                                            window,
                                            CAlignmentSmear::eSmearStrand_Both,
                                            seq_annot));
            if (ref) {
                aligns.push_back(ref);
            }
        }
    }
}


static bool
s_CheckTraceMateAlign(const CSeq_align& align,
                      int& ti,
                      int& mate_ti,
                      int expect_mate_ti = 0)
{
    //
    // specific requirements for mate pair alignments:
    //

    // must be a pairwise dense-seg
    if ( !align.GetSegs().IsDenseg()  ||
         align.GetSegs().GetDenseg().GetIds().size() != 2) {
        return false;
    }

    // check the IDs for a trace ID
    ti = 0;
    ITERATE (CDense_seg::TIds, iter,
             align.GetSegs().GetDenseg().GetIds()) {
        const CSeq_id& id = **iter;
        if (id.IsGeneral()  &&
            (id.GetGeneral().GetDb() == "ti"  ||
            id.GetGeneral().GetDb() == "TRACE")) {
            ti = id.GetGeneral().GetTag().GetId();
            break;
        }
    }

    if (ti == 0) {
        return false;
    }

    // must have a score field named 'matepair_ti'
    mate_ti = 0;
    if ( !align.GetNamedScore("matepair ti", mate_ti)) {
        if ( !align.GetNamedScore("bad matepair ti", mate_ti)) {
            return false;
        }
    }

    if (expect_mate_ti  &&  mate_ti != expect_mate_ti) {
        return false;
    }

    return true;
}


//
// GetAlignments()
// this retrieves a set of alignments from the document
//
void CSeqUtils::GetAlignments(const CBioseq_Handle&    handle,
                              const TSeqRange&         range,
                              CLayoutEngine::TObjects& aligns,
                              TAlignFlags flags)
{
    aligns.clear();
    if ( !handle ) {
        return;
    }

    SAnnotSelector selector = GetAnnotSelector(CSeq_annot::TData::e_Align);
    CAlign_CI align_iter(handle, range.GetFrom(), range.GetTo(), selector);

    typedef pair<CRef<CLayoutPWAlign> , int> TMatedAlign;
    typedef map<int, TMatedAlign> TMatedAlignments;
    TMatedAlignments mated_aligns;

    for (;  align_iter ;  ++align_iter) {
        const CSeq_align& align = *align_iter;

        try {

            // generate an alignment manager for this alignment
            CRef<CAlnVec> aln_mgr;
            if (align.GetSegs().IsDenseg()) {
                aln_mgr.Reset(new CAlnVec(align.GetSegs().GetDenseg(),
                                          handle.GetScope()));
            } else {
                CAlnMix mix(handle.GetScope());
                mix.Add(align);
                mix.Merge(CAlnMix::fGen2EST |
                          CAlnMix::fTryOtherMethodOnFail |
                          CAlnMix::fGapJoin);

                aln_mgr.Reset(new CAlnVec(mix.GetDenseg(), handle.GetScope()));
            }

            // anchor on the referent sequence
            CAlnVec::TNumrow row = 0;
            CAlnVec::TNumrow anchor = 0;
            for (row = 0;  row != aln_mgr->GetNumRows();  ++row) {
                if ( handle.IsSynonym(aln_mgr->GetSeqId(row)) ) {
                    aln_mgr->SetAnchor(row);
                    anchor = row;
                    break;
                }
            }

            //
            // check to see if this is a mate pair alignment
            //
            int ti = 0;
            int mate_ti = 0;
            if ( s_CheckTraceMateAlign(align, ti, mate_ti) ) {
                CRef<CLayoutPWAlign> pwal
                    (new CLayoutPWAlign(*aln_mgr, align));
                if (flags & fAlign_LinkMatePairs) {
                    // we will process mate mair alignments in a separate pass
                    // after we complete all other alignments
                    TMatedAlign mp(pwal, mate_ti);
                    mated_aligns[ti] = mp;
                } else {
                    CLayoutMatePair::TAlignList als;
                    als.push_back(pwal);
                    CRef<CLayoutObject> obj(new CLayoutMatePair(als));
                    aligns.push_back(obj);
                }
            } else {
                // "normal" alignment
                CRef<CLayoutObject> ref;
                if (aln_mgr->GetNumRows() == 2) {
                    ref.Reset(new CLayoutPWAlign(*aln_mgr, align));
                } else {
                    ref.Reset(new CLayoutAlign(*aln_mgr, align));
                }
                aligns.push_back(ref);
            }
        }
        catch (CException& e) {
            // errors ignored
            _TRACE("error in GetAlignments(): " << e.GetMsg());
        }
    }

    //
    // final pass - connect our mate pair alignments, if we can
    //
    NON_CONST_ITERATE (TMatedAlignments, iter, mated_aligns) {
        if ( !iter->second.first ) {
            continue;
        }

        CLayoutPWAlign& first_mate    = *iter->second.first;
        const CSeq_align&     align   = first_mate.GetAlignment();
        int                   ti      = iter->first;
        int                   mate_ti = iter->second.second;

        //
        // verify that we have a mated alignment for this trace alignment
        //
        CRef<CLayoutPWAlign> second_mate;
        TMatedAlignments::iterator mate_iter = mated_aligns.find(mate_ti);
        if (mate_iter == mated_aligns.end()) {
            //
            // NOT FOUND
            // try iterating
            //
            CSeq_id id("gnl|ti|" + NStr::IntToString(mate_ti));
            CBioseq_Handle h = 
                handle.GetScope().GetBioseqHandle(id);
            CAlign_CI second_mate_iter(h, 0, h.GetBioseqLength(),
                                       selector);
            for ( ;  second_mate_iter;  ++second_mate_iter) {
                int temp = 0;
                if ( !s_CheckTraceMateAlign(*second_mate_iter, temp,
                                            mate_ti, ti) ) {
                    continue;
                }

                CRef<CAlnVec> aln_mgr(new CAlnVec(align.GetSegs().GetDenseg(),
                                                  handle.GetScope()));

                // anchor on the referent sequence
                CAlnVec::TNumrow row = 0;
                CAlnVec::TNumrow anchor = 0;
                for (row = 0;  row != aln_mgr->GetNumRows();  ++row) {
                    if ( handle.IsSynonym(aln_mgr->GetSeqId(row)) ) {
                        aln_mgr->SetAnchor(row);
                        anchor = row;
                        break;
                    }
                }

                second_mate.Reset(new CLayoutPWAlign(*aln_mgr, align));
                break;
            }

            if ( !second_mate ) {
                LOG_POST(Error << "failed to find mate pair for ti = " << ti);

                // no mate available, so treat it as a normal alignment
                // we need to trim leading and trailing gaps to get a good
                // placement
                CRef<CLayoutObject> ref
                    (new CLayoutPWAlign(first_mate.GetAlignMgr(), align));
                aligns.push_back(ref);
                continue;
            }
        } else if ( !mate_iter->second.first ) {
            // already processed
            continue;
        } else {
            second_mate = mate_iter->second.first;
            mate_iter->second.first.Reset();
        }