align_to_neighbors.cpp

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

/*
 * ===========================================================================
 * PRODUCTION $Log: align_to_neighbors.cpp,v $
 * PRODUCTION Revision 1000.4  2004/06/01 20:54:09  gouriano
 * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.17
 * PRODUCTION
 * ===========================================================================
 */

/*  $Id: align_to_neighbors.cpp,v 1000.4 2004/06/01 20:54:09 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:  Josh Cherry
 *
 * File Description:
 *    gbench plugin for aligning to neighbors
 */

#include <ncbi_pch.hpp>
#include "align_to_neighbors.hpp"
#include "blast_util.hpp"

#include <objects/entrez2/entrez2_client.hpp>

#include <gui/core/idocument.hpp>
#include <gui/core/plugin_utils.hpp>
#include <gui/plugin/PluginCommandSet.hpp>
#include <gui/plugin/PluginValue.hpp>
#include <gui/plugin/PluginInfo.hpp>
#include <gui/plugin/PluginRequest.hpp>
#include <gui/core/version.hpp>
#include <objmgr/util/feature.hpp>
#include <objmgr/util/sequence.hpp>
#include <gui/objutils/utils.hpp>
#include <gui/core/doc_manager.hpp>
#include <gui/utils/message_box.hpp>

#include <objmgr/scope.hpp>
#include <objmgr/seq_vector.hpp>

#include <objects/seqloc/Seq_loc.hpp>
#include <objects/seqloc/Seq_point.hpp>
#include <objects/seqfeat/Gb_qual.hpp>

#include <algo/blast/api/bl2seq.hpp>

BEGIN_NCBI_SCOPE
USING_SCOPE(ncbi::objects);


void CAlignToNeighbors::Run(CPluginMessage& msg)
{
    CProgressDlgEx progress;
    progress.Show();
    progress.SetTitle("Computing Alignment to Entrez Neighbors");

    const CPluginCommand& cmd  = msg.GetRequest().GetCommand();
    CPluginReply&         reply = msg.SetReply();
    reply.SetStatus(eMessageStatus_failed);

    const CSeq_loc*  loc = NULL;
    const CSeq_loc*  parent_loc = NULL;
    const IDocument* doc = NULL;

    blast::EProgram prog = blast::eBlastn;
    string db;

    if (cmd.HasArgument("mrna_feat")  &&
        CPluginUtils::IsValid(cmd["mrna_feat"])) {

        const CPluginArg& arg = cmd["mrna_feat"];
        doc = arg.GetDocument();
        const CSeq_feat*  feat =
            dynamic_cast<const CSeq_feat*>(arg.GetObject());
        loc = &feat->GetProduct();

        prog = blast::eBlastn;

        if (cmd["genomic"].AsBoolean()) {
            parent_loc = &feat->GetLocation();
        }
        db = "nucleotide";

    } else if (cmd.HasArgument("cds_feat")  &&
               CPluginUtils::IsValid(cmd["cds_feat"])) {

        const CPluginArg& arg = cmd["cds_feat"];
        doc = arg.GetDocument();
        const CSeq_feat*  feat =
            dynamic_cast<const CSeq_feat*>(arg.GetObject());
        loc = &feat->GetProduct();

        prog = blast::eBlastp;
        db = "protein";
    }

    if ( !loc  ||  !doc ) {
        NcbiMessageBox("Invalid inputs");
        return;
    }

    CScope& scope = doc->GetScope();

    //
    // retrieve the 'gi' of the sequence
    //
    int gi = 0;
    {{
        const CSeq_id& id = sequence::GetId(*loc);
        if (id.IsGi()) {
            gi = id.GetGi();
        }
    }}

    if (gi == 0) {
        try {
            CBioseq_Handle handle = scope.GetBioseqHandle(*loc);
            if (handle) {
                const CSeq_id& gi_id =
                    sequence::GetId(handle, sequence::eGetId_ForceGi);
                gi = gi_id.GetGi();
            }
        }
        catch (...) {
            // not a GI id
        }
    }

    if (gi == 0) {
        NcbiMessageBox("The selected sequence doesn't have a valid "
                       " GenBank ID");
        return;
    }

    progress.SetMessage("Retrieving Entrez neighbors...");

    // Get ids of nucleotide neighbors
    vector<int> neighbor_gis;
    CEntrez2Client e2c;
    e2c.GetNeighbors(gi, db, db, neighbor_gis);

    if (neighbor_gis.empty()) {
        NcbiMessageBox(string("No neighbors found for gi|") +
                       NStr::IntToString(gi));
        return;
    }

    progress.SetMessage("Found " + NStr::IntToString(neighbor_gis.size()) +
                      " neighbors...");
    _TRACE("Found " << neighbor_gis.size() << " neighbors");

    // Neighbor filtering

    // First defend against an easy mistake to make
    if (cmd.HasArgument("query_string")
        && CPluginUtils::IsValid(cmd["query_string"])) {
        if (!(cmd.HasArgument("neighbor_filter")
              && CPluginUtils::IsValid(cmd["neighbor_filter"])
              && cmd["neighbor_filter"].AsString() ==
              "arbitrary query string (below)")) {
            NcbiMessageBox("Query string filled in, "
                           "but that type of filtering not selected");
            return;
        }
    }

    // Now filter neighbors if requested by user
    if (cmd.HasArgument("neighbor_filter")
        && CPluginUtils::IsValid(cmd["neighbor_filter"])) {
        string filter = cmd["neighbor_filter"].AsString();
        if (filter != "all") {
            string query_string;
            if (filter == "mRNA only") {
                query_string = "biomol_mrna[Prop]";
            } else if (filter == "genomic only") {
                query_string = "biomol_genomic[Prop]";
            } else if (filter == "arbitrary query string (below)") {
                if (!(cmd.HasArgument("query_string")
                      && CPluginUtils::IsValid(cmd["query_string"]))) {
                    NcbiMessageBox("Filter by string selected, "
                                   "but query string not filled in");
                    return;
                }
                query_string = cmd["query_string"].AsString();
            } else {
                // this shouldn't happen
            }

            progress.SetMessage("Filtering neighbors...");

            vector<int> filtered_gis;
            e2c.FilterIds(neighbor_gis, db, query_string, filtered_gis);
            neighbor_gis.swap(filtered_gis);
        }
    }

    if (neighbor_gis.empty()) {
        NcbiMessageBox("No neighbors matched filter conditions");
        return;
    }

    _TRACE(neighbor_gis.size() << " neighbors matched filtering conditions");

    //
    // retrieve our sequences
    // technically, this is not necessary, but we do it for GUI appeal...
    //
    {{
        int counter = 0;
        ITERATE (vector<int>, iter, neighbor_gis) {
            CSeq_id id;
            id.SetGi(*iter);
            CBioseq_Handle handle = scope.GetBioseqHandle(id);
            const CSeq_id& best_id = sequence::GetId(handle,
                                                     sequence::eGetId_Best);

            ++counter;
            string msg = "Retrieved ";
            best_id.GetLabel(&msg);
            msg += " (" +
                NStr::IntToString(counter) +
                "/" + NStr::IntToString(neighbor_gis.size()) + "), " +
                NStr::IntToString(handle.GetBioseqLength()) + " bases";
            progress.SetMessage(msg);
            progress.SetPctCompleted(100 * counter / neighbor_gis.size());

            /**
            if (x_IsInterrupted()) {
                return;
            }
            **/
        }
    }}

    // Construct CBl2Seq object

    progress.SetMessage("Preparing to BLAST...");

    blast::SSeqLoc bl_query(loc, &scope);
    blast::TSeqLocVector targets;
    if (parent_loc) {
        // we blast everything in the context of the parent
        targets.push_back (bl_query);
        bl_query = blast::SSeqLoc(parent_loc, &scope);
    }

    ITERATE (vector<int>, iter, neighbor_gis) {
        CRef<CSeq_loc> subject_loc(new CSeq_loc());
        subject_loc->SetWhole().SetGi(*iter);
        blast::SSeqLoc bl_subject(subject_loc, &scope);
        targets.push_back(bl_subject);
    }

    blast::CBl2Seq blaster(bl_query, targets, prog);
    CBlastUtils::ArgsToBlastOptions(cmd, blaster.SetOptions());

    // run blast
    progress.SetMessage("Computing alignment...");
	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("Neighbor alignment");

    /**
    // final check for completion
    if (x_IsInterrupted()) {
        return;
    }
    **/

    // standard reply handling
    reply.AddObject(*doc, *annot);
    reply.AddAction(CPluginReplyAction::e_Add_to_document);
    reply.SetStatus(eMessageStatus_success);

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

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


END_NCBI_SCOPE


/*
 * ===========================================================================
 * $Log: align_to_neighbors.cpp,v $
 * Revision 1000.4  2004/06/01 20:54:09  gouriano
 * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.17
 *
 * Revision 1.17  2004/05/21 22:27:46  gorelenk
 * Added PCH ncbi_pch.hpp
 *
 * Revision 1.16  2004/05/20 12:32:45  dicuccio
 * Use GetBioseqLength() instead of CSeqVector::Size()
 *
 * Revision 1.15  2004/05/15 03:17:04  ucko
 * Add missing #includes (formerly indirect?)
 *
 * Revision 1.14  2004/05/03 13:05:42  dicuccio
 * gui/utils --> gui/objutils where needed
 *
 * Revision 1.13  2004/03/05 17:33:44  dicuccio
 * Changed plugin to run single-threaded until thread issues can be sorted out.
 * Use sequence::GetId() instead of CSeq_id::GetStringDescr()
 *
 * Revision 1.12  2004/01/27 18:40:02  dicuccio
 * Code clean-up.  Renamed plugin classes to follow standard pattern
 *
 * Revision 1.11  2004/01/13 20:36:29  dicuccio
 * Use CBlastUtils for standard argument processing.  Make sure to pass the new
 * objects back to the framework with the appropriate action set (add to the
 * document specified)
 *
 * Revision 1.10  2003/12/04 18:13:23  dicuccio
 * Changed to match API change in CProgressDlgEx
 *
 * Revision 1.9  2003/11/26 21:08:07  ucko
 * Adjust for current CBlastOptions location and API.
 *
 * Revision 1.8  2003/11/26 17:12:30  dicuccio
 * Switched to use CThreadedAlgorithm<>
 *
 * Revision 1.7  2003/11/04 17:49:22  dicuccio
 * Changed calling parameters for plugins - pass CPluginMessage instead of paired
 * CPluginCommand/CPluginReply
 *
 * Revision 1.6  2003/11/03 17:41:19  dicuccio
 * Fixed to match changes in BLAST API
 *
 * Revision 1.5  2003/10/23 16:21:28  dicuccio
 * Remvoed dead code
 *
 * Revision 1.4  2003/10/16 21:52:07  jcherry
 * Added filtering of neighbors
 *
 * ===========================================================================
 */