align_5prime.cpp

来自「ncbi源码」· C++ 代码 · 共 252 行

/*
 * ===========================================================================
 * PRODUCTION $Log: align_5prime.cpp,v $
 * PRODUCTION Revision 1000.1  2004/06/01 20:54:04  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: align_5prime.cpp,v 1000.1 2004/06/01 20:54:04 gouriano Exp $
 * ===========================================================================
 *
 *                            PUBLIC DOMAIN NOTICE
 *               National Center for Biotechnology Information
 *
 *  This software/database is a "United States Government Work" under the
 *  terms of the United States Copyright Act.  It was written as part of
 *  the author's official duties as a United States Government employee and
 *  thus cannot be copyrighted.  This software/database is freely available
 *  to the public for use. The National Library of Medicine and the U.S.
 *  Government have not placed any restriction on its use or reproduction.
 *
 *  Although all reasonable efforts have been taken to ensure the accuracy
 *  and reliability of the software and data, the NLM and the U.S.
 *  Government do not and cannot warrant the performance or results that
 *  may be obtained by using this software or data. The NLM and the U.S.
 *  Government disclaim all warranties, express or implied, including
 *  warranties of performance, merchantability or fitness for any particular
 *  purpose.
 *
 *  Please cite the author in any work or product based on this material.
 *
 * ===========================================================================
 *
 * Authors:  Mike DiCuccio
 *
 * File Description:
 *    gbench plugin for aligning to neighbors
 */

#include <ncbi_pch.hpp>
#include "align_5prime.hpp"

#include "blast_util.hpp"
#include <algo/blast/api/bl2seq.hpp>
#include <gui/core/doc_manager.hpp>
#include <gui/core/plugin_utils.hpp>
#include <gui/core/version.hpp>
#include <gui/plugin/AlgoCommand.hpp>
#include <gui/plugin/PluginCommand.hpp>
#include <gui/plugin/PluginCommandSet.hpp>
#include <gui/plugin/PluginInfo.hpp>
#include <gui/plugin/PluginReply.hpp>
#include <gui/plugin/PluginRequest.hpp>
#include <gui/plugin/PluginValueConstraint.hpp>
#include <objects/seqalign/Seq_align.hpp>
#include <objects/seqfeat/Seq_feat.hpp>
#include <objects/seqloc/Seq_loc.hpp>
#include <objmgr/util/feature.hpp>
#include <objmgr/util/sequence.hpp>


BEGIN_NCBI_SCOPE
USING_SCOPE(ncbi::objects);


void CAlgoPlugin_Align5Prime::GetInfo(CPluginInfo& info)
{
    info.Reset();

    // version info macro
    info.SetInfo(CPluginVersion::eMajor, CPluginVersion::eMinor, 0,
                 string(__DATE__) + " " + string(__TIME__),
                 "CAlgoPlugin_Align5Prime",
                 "Alignments/Align 5' region of genes",
                 "Create an alignment of the 5' flanking regions of genes",
                 "");

    // command info
    CPluginCommandSet& cmds = info.SetCommands();
    CPluginCommand&    args = cmds.AddAlgoCommand(eAlgoCommand_run);

    args.AddArgument("genes", "Genes to align",
                     CSeq_feat::GetTypeInfo(),
                     CPluginArg::TData::e_Array);
    args.SetConstraint("genes",
                       (*CPluginValueConstraint::CreateFeatSubtype(),
                        CSeqFeatData::eSubtype_gene));
    args.AddDefaultArgument("offset", "Alignment offset",
                            CPluginArg::eInteger, "10000");
    args.AddDefaultArgument("length", "Alignment length",
                            CPluginArg::eInteger, "10000");

    CBlastUtils::AddBlastArgs(args, blast::eBlastn);
}


void CAlgoPlugin_Align5Prime::RunCommand(CPluginMessage& msg)
{
    const CPluginCommand& args = msg.GetRequest().GetCommand();
    CPluginReply&         reply = msg.SetReply();
    reply.SetStatus(eMessageStatus_failed);

    plugin_args::TFeatList feats;
    GetArgValue(args["genes"], feats);
    TSeqPos offs = args["offset"].AsInteger();
    TSeqPos len  = args["offset"].AsInteger();

    // scan our features.  We must make sure that we get the appropriate
    // upstream location - some features may be too long
    blast::TSeqLocVector targets;
    ITERATE (plugin_args::TFeatList, iter, feats) {
        const CSeq_feat& feat = *iter->second;
        const IDocument& doc  = *iter->first;

        CRef<CSeq_loc> loc(new CSeq_loc());
        ENa_strand strand = sequence::GetStrand(feat.GetLocation());
        TSeqRange range = feat.GetLocation().GetTotalRange();

        TSeqPos from;
        TSeqPos to;
        if (strand == eNa_strand_minus) {
            from = range.GetFrom() + offs;

            // make sure we're bounded by the sequence upstream
            CBioseq_Handle handle =
                doc.GetScope().GetBioseqHandle(feat.GetLocation());
            if (from > handle.GetBioseqLength()) {
                from = handle.GetBioseqLength();
            }

            to = from - len;
            if (len > from) {
                to = 0;
            }

        } else {
            // assume '+' strand
            from = range.GetFrom() - offs;

            // make sure we're bounded by the sequence upstream
            if (range.GetFrom() < offs) {
                from = 0;
            }

            to = from + len;

            CBioseq_Handle handle =
                doc.GetScope().GetBioseqHandle(feat.GetLocation());
            if (to > handle.GetBioseqLength()) {
                to = handle.GetBioseqLength();
            }
        }

        loc->SetInt().SetFrom(from);
        loc->SetInt().SetTo  (to);
        loc->SetInt().SetStrand(strand);
        loc->SetId(sequence::GetId(feat.GetLocation()));
        blast::SSeqLoc blast_loc(loc, &doc.GetScope());
        targets.push_back(blast_loc);
    }

    // find a long one to be the query
    blast::TSeqLocVector::iterator iter = targets.begin();
    blast::TSeqLocVector::iterator best = targets.begin();
    ++iter;
    for ( ;  iter != targets.end();  ++iter) {
        if (iter->seqloc->GetTotalRange().GetLength() >
            best->seqloc->GetTotalRange().GetLength()) {
            best = iter;
        }
    }

    blast::SSeqLoc query(*best);
    targets.erase(best);

    blast::CBl2Seq blaster(query, targets, blast::eBlastn);
    CBlastUtils::ArgsToBlastOptions(args, blaster.SetOptions());

    // run blast
	blast::TSeqAlignVector aligns = blaster.Run();

    // make an annotation
    CRef<CSeq_annot> annot(new CSeq_annot());
	ITERATE (blast::TSeqAlignVector, iter, aligns) {
		annot->SetData().SetAlign()
			.insert(annot->SetData().SetAlign().end(),
					(*iter)->Get().begin(), (*iter)->Get().end());
	}
    annot->AddName("Upstream flanking alignment");

    //
    // pass back to the system.  We may use the same scope and just attach,
    // if that is appropriate
    //
    CConstRef<IDocument> doc_ref(feats.begin()->first);
    ITERATE (plugin_args::TFeatList, iter, feats) {
        if ( !doc_ref ) {
            doc_ref.Reset(iter->first);
        } else if (iter->first != doc_ref) {
            doc_ref.Reset();
            break;
        }
    }

    if ( !doc_ref ) {
        //
        // features come from different documents
        // create a new one to handle the results
        //
        CRef<CScope> new_scope(new CScope(CDocManager::GetObjectManager()));
        ITERATE (plugin_args::TFeatList, iter, feats) {
            new_scope->AddScope(iter->first->GetScope());
        }

        doc_ref.Reset(CDocManager::CreateDocument(*new_scope, *annot));
    } else {
        reply.AddAction(CPluginReplyAction::e_Add_to_document);
    }

    // standard reply handling
    reply.AddObject(*doc_ref, *annot);
    reply.SetStatus(eMessageStatus_success);

    // go ahead and create an alignment view
    CRef<CSelectionBuffer> buf(new CSelectionBuffer());
    buf->AddSelection(doc_ref, *annot);

    CPluginUtils::CallPlugin("CAlnMultiView",
                             eViewCommand_new_view, *buf);
}


END_NCBI_SCOPE


/*
 * ===========================================================================
 * $Log: align_5prime.cpp,v $
 * Revision 1000.1  2004/06/01 20:54:04  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
 *
 * Revision 1.2  2004/05/21 22:27:46  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.1  2004/01/27 18:39:01  dicuccio
 * Initial revision
 *
 * Revision 1.11  2004/01/13 20:36:29  dicuccio
 * ===========================================================================
 */