align_to_neighbors.cpp
上传用户:yhdzpy8989
上传日期:2007-06-13
资源大小:13604k
文件大小:12k
源码类别:

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: align_to_neighbors.cpp,v $
  4.  * PRODUCTION Revision 1000.4  2004/06/01 20:54:09  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.17
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: align_to_neighbors.cpp,v 1000.4 2004/06/01 20:54:09 gouriano Exp $
  10.  * ===========================================================================
  11.  *
  12.  *                            PUBLIC DOMAIN NOTICE
  13.  *               National Center for Biotechnology Information
  14.  *
  15.  *  This software/database is a "United States Government Work" under the
  16.  *  terms of the United States Copyright Act.  It was written as part of
  17.  *  the author's official duties as a United States Government employee and
  18.  *  thus cannot be copyrighted.  This software/database is freely available
  19.  *  to the public for use. The National Library of Medicine and the U.S.
  20.  *  Government have not placed any restriction on its use or reproduction.
  21.  *
  22.  *  Although all reasonable efforts have been taken to ensure the accuracy
  23.  *  and reliability of the software and data, the NLM and the U.S.
  24.  *  Government do not and cannot warrant the performance or results that
  25.  *  may be obtained by using this software or data. The NLM and the U.S.
  26.  *  Government disclaim all warranties, express or implied, including
  27.  *  warranties of performance, merchantability or fitness for any particular
  28.  *  purpose.
  29.  *
  30.  *  Please cite the author in any work or product based on this material.
  31.  *
  32.  * ===========================================================================
  33.  *
  34.  * Authors:  Josh Cherry
  35.  *
  36.  * File Description:
  37.  *    gbench plugin for aligning to neighbors
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include "align_to_neighbors.hpp"
  41. #include "blast_util.hpp"
  42. #include <objects/entrez2/entrez2_client.hpp>
  43. #include <gui/core/idocument.hpp>
  44. #include <gui/core/plugin_utils.hpp>
  45. #include <gui/plugin/PluginCommandSet.hpp>
  46. #include <gui/plugin/PluginValue.hpp>
  47. #include <gui/plugin/PluginInfo.hpp>
  48. #include <gui/plugin/PluginRequest.hpp>
  49. #include <gui/core/version.hpp>
  50. #include <objmgr/util/feature.hpp>
  51. #include <objmgr/util/sequence.hpp>
  52. #include <gui/objutils/utils.hpp>
  53. #include <gui/core/doc_manager.hpp>
  54. #include <gui/utils/message_box.hpp>
  55. #include <objmgr/scope.hpp>
  56. #include <objmgr/seq_vector.hpp>
  57. #include <objects/seqloc/Seq_loc.hpp>
  58. #include <objects/seqloc/Seq_point.hpp>
  59. #include <objects/seqfeat/Gb_qual.hpp>
  60. #include <algo/blast/api/bl2seq.hpp>
  61. BEGIN_NCBI_SCOPE
  62. USING_SCOPE(ncbi::objects);
  63. void CAlignToNeighbors::Run(CPluginMessage& msg)
  64. {
  65.     CProgressDlgEx progress;
  66.     progress.Show();
  67.     progress.SetTitle("Computing Alignment to Entrez Neighbors");
  68.     const CPluginCommand& cmd  = msg.GetRequest().GetCommand();
  69.     CPluginReply&         reply = msg.SetReply();
  70.     reply.SetStatus(eMessageStatus_failed);
  71.     const CSeq_loc*  loc = NULL;
  72.     const CSeq_loc*  parent_loc = NULL;
  73.     const IDocument* doc = NULL;
  74.     blast::EProgram prog = blast::eBlastn;
  75.     string db;
  76.     if (cmd.HasArgument("mrna_feat")  &&
  77.         CPluginUtils::IsValid(cmd["mrna_feat"])) {
  78.         const CPluginArg& arg = cmd["mrna_feat"];
  79.         doc = arg.GetDocument();
  80.         const CSeq_feat*  feat =
  81.             dynamic_cast<const CSeq_feat*>(arg.GetObject());
  82.         loc = &feat->GetProduct();
  83.         prog = blast::eBlastn;
  84.         if (cmd["genomic"].AsBoolean()) {
  85.             parent_loc = &feat->GetLocation();
  86.         }
  87.         db = "nucleotide";
  88.     } else if (cmd.HasArgument("cds_feat")  &&
  89.                CPluginUtils::IsValid(cmd["cds_feat"])) {
  90.         const CPluginArg& arg = cmd["cds_feat"];
  91.         doc = arg.GetDocument();
  92.         const CSeq_feat*  feat =
  93.             dynamic_cast<const CSeq_feat*>(arg.GetObject());
  94.         loc = &feat->GetProduct();
  95.         prog = blast::eBlastp;
  96.         db = "protein";
  97.     }
  98.     if ( !loc  ||  !doc ) {
  99.         NcbiMessageBox("Invalid inputs");
  100.         return;
  101.     }
  102.     CScope& scope = doc->GetScope();
  103.     //
  104.     // retrieve the 'gi' of the sequence
  105.     //
  106.     int gi = 0;
  107.     {{
  108.         const CSeq_id& id = sequence::GetId(*loc);
  109.         if (id.IsGi()) {
  110.             gi = id.GetGi();
  111.         }
  112.     }}
  113.     if (gi == 0) {
  114.         try {
  115.             CBioseq_Handle handle = scope.GetBioseqHandle(*loc);
  116.             if (handle) {
  117.                 const CSeq_id& gi_id =
  118.                     sequence::GetId(handle, sequence::eGetId_ForceGi);
  119.                 gi = gi_id.GetGi();
  120.             }
  121.         }
  122.         catch (...) {
  123.             // not a GI id
  124.         }
  125.     }
  126.     if (gi == 0) {
  127.         NcbiMessageBox("The selected sequence doesn't have a valid "
  128.                        " GenBank ID");
  129.         return;
  130.     }
  131.     progress.SetMessage("Retrieving Entrez neighbors...");
  132.     // Get ids of nucleotide neighbors
  133.     vector<int> neighbor_gis;
  134.     CEntrez2Client e2c;
  135.     e2c.GetNeighbors(gi, db, db, neighbor_gis);
  136.     if (neighbor_gis.empty()) {
  137.         NcbiMessageBox(string("No neighbors found for gi|") +
  138.                        NStr::IntToString(gi));
  139.         return;
  140.     }
  141.     progress.SetMessage("Found " + NStr::IntToString(neighbor_gis.size()) +
  142.                       " neighbors...");
  143.     _TRACE("Found " << neighbor_gis.size() << " neighbors");
  144.     // Neighbor filtering
  145.     // First defend against an easy mistake to make
  146.     if (cmd.HasArgument("query_string")
  147.         && CPluginUtils::IsValid(cmd["query_string"])) {
  148.         if (!(cmd.HasArgument("neighbor_filter")
  149.               && CPluginUtils::IsValid(cmd["neighbor_filter"])
  150.               && cmd["neighbor_filter"].AsString() ==
  151.               "arbitrary query string (below)")) {
  152.             NcbiMessageBox("Query string filled in, "
  153.                            "but that type of filtering not selected");
  154.             return;
  155.         }
  156.     }
  157.     // Now filter neighbors if requested by user
  158.     if (cmd.HasArgument("neighbor_filter")
  159.         && CPluginUtils::IsValid(cmd["neighbor_filter"])) {
  160.         string filter = cmd["neighbor_filter"].AsString();
  161.         if (filter != "all") {
  162.             string query_string;
  163.             if (filter == "mRNA only") {
  164.                 query_string = "biomol_mrna[Prop]";
  165.             } else if (filter == "genomic only") {
  166.                 query_string = "biomol_genomic[Prop]";
  167.             } else if (filter == "arbitrary query string (below)") {
  168.                 if (!(cmd.HasArgument("query_string")
  169.                       && CPluginUtils::IsValid(cmd["query_string"]))) {
  170.                     NcbiMessageBox("Filter by string selected, "
  171.                                    "but query string not filled in");
  172.                     return;
  173.                 }
  174.                 query_string = cmd["query_string"].AsString();
  175.             } else {
  176.                 // this shouldn't happen
  177.             }
  178.             progress.SetMessage("Filtering neighbors...");
  179.             vector<int> filtered_gis;
  180.             e2c.FilterIds(neighbor_gis, db, query_string, filtered_gis);
  181.             neighbor_gis.swap(filtered_gis);
  182.         }
  183.     }
  184.     if (neighbor_gis.empty()) {
  185.         NcbiMessageBox("No neighbors matched filter conditions");
  186.         return;
  187.     }
  188.     _TRACE(neighbor_gis.size() << " neighbors matched filtering conditions");
  189.     //
  190.     // retrieve our sequences
  191.     // technically, this is not necessary, but we do it for GUI appeal...
  192.     //
  193.     {{
  194.         int counter = 0;
  195.         ITERATE (vector<int>, iter, neighbor_gis) {
  196.             CSeq_id id;
  197.             id.SetGi(*iter);
  198.             CBioseq_Handle handle = scope.GetBioseqHandle(id);
  199.             const CSeq_id& best_id = sequence::GetId(handle,
  200.                                                      sequence::eGetId_Best);
  201.             ++counter;
  202.             string msg = "Retrieved ";
  203.             best_id.GetLabel(&msg);
  204.             msg += " (" +
  205.                 NStr::IntToString(counter) +
  206.                 "/" + NStr::IntToString(neighbor_gis.size()) + "), " +
  207.                 NStr::IntToString(handle.GetBioseqLength()) + " bases";
  208.             progress.SetMessage(msg);
  209.             progress.SetPctCompleted(100 * counter / neighbor_gis.size());
  210.             /**
  211.             if (x_IsInterrupted()) {
  212.                 return;
  213.             }
  214.             **/
  215.         }
  216.     }}
  217.     // Construct CBl2Seq object
  218.     progress.SetMessage("Preparing to BLAST...");
  219.     blast::SSeqLoc bl_query(loc, &scope);
  220.     blast::TSeqLocVector targets;
  221.     if (parent_loc) {
  222.         // we blast everything in the context of the parent
  223.         targets.push_back (bl_query);
  224.         bl_query = blast::SSeqLoc(parent_loc, &scope);
  225.     }
  226.     ITERATE (vector<int>, iter, neighbor_gis) {
  227.         CRef<CSeq_loc> subject_loc(new CSeq_loc());
  228.         subject_loc->SetWhole().SetGi(*iter);
  229.         blast::SSeqLoc bl_subject(subject_loc, &scope);
  230.         targets.push_back(bl_subject);
  231.     }
  232.     blast::CBl2Seq blaster(bl_query, targets, prog);
  233.     CBlastUtils::ArgsToBlastOptions(cmd, blaster.SetOptions());
  234.     // run blast
  235.     progress.SetMessage("Computing alignment...");
  236. blast::TSeqAlignVector aligns = blaster.Run();
  237.     // make an annotation
  238.     CRef<CSeq_annot> annot(new CSeq_annot());
  239. ITERATE (blast::TSeqAlignVector, iter, aligns) {
  240. annot->SetData().SetAlign()
  241. .insert(annot->SetData().SetAlign().end(),
  242. (*iter)->Get().begin(), (*iter)->Get().end());
  243. }
  244.     annot->AddName("Neighbor alignment");
  245.     /**
  246.     // final check for completion
  247.     if (x_IsInterrupted()) {
  248.         return;
  249.     }
  250.     **/
  251.     // standard reply handling
  252.     reply.AddObject(*doc, *annot);
  253.     reply.AddAction(CPluginReplyAction::e_Add_to_document);
  254.     reply.SetStatus(eMessageStatus_success);
  255.     // go ahead and create an alignment view
  256.     CRef<CSelectionBuffer> buf(new CSelectionBuffer());
  257.     buf->AddSelection(doc, *annot);
  258.     CPluginUtils::CallPlugin("CAlnMultiView",
  259.                              eViewCommand_new_view, *buf);
  260. }
  261. END_NCBI_SCOPE
  262. /*
  263.  * ===========================================================================
  264.  * $Log: align_to_neighbors.cpp,v $
  265.  * Revision 1000.4  2004/06/01 20:54:09  gouriano
  266.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.17
  267.  *
  268.  * Revision 1.17  2004/05/21 22:27:46  gorelenk
  269.  * Added PCH ncbi_pch.hpp
  270.  *
  271.  * Revision 1.16  2004/05/20 12:32:45  dicuccio
  272.  * Use GetBioseqLength() instead of CSeqVector::Size()
  273.  *
  274.  * Revision 1.15  2004/05/15 03:17:04  ucko
  275.  * Add missing #includes (formerly indirect?)
  276.  *
  277.  * Revision 1.14  2004/05/03 13:05:42  dicuccio
  278.  * gui/utils --> gui/objutils where needed
  279.  *
  280.  * Revision 1.13  2004/03/05 17:33:44  dicuccio
  281.  * Changed plugin to run single-threaded until thread issues can be sorted out.
  282.  * Use sequence::GetId() instead of CSeq_id::GetStringDescr()
  283.  *
  284.  * Revision 1.12  2004/01/27 18:40:02  dicuccio
  285.  * Code clean-up.  Renamed plugin classes to follow standard pattern
  286.  *
  287.  * Revision 1.11  2004/01/13 20:36:29  dicuccio
  288.  * Use CBlastUtils for standard argument processing.  Make sure to pass the new
  289.  * objects back to the framework with the appropriate action set (add to the
  290.  * document specified)
  291.  *
  292.  * Revision 1.10  2003/12/04 18:13:23  dicuccio
  293.  * Changed to match API change in CProgressDlgEx
  294.  *
  295.  * Revision 1.9  2003/11/26 21:08:07  ucko
  296.  * Adjust for current CBlastOptions location and API.
  297.  *
  298.  * Revision 1.8  2003/11/26 17:12:30  dicuccio
  299.  * Switched to use CThreadedAlgorithm<>
  300.  *
  301.  * Revision 1.7  2003/11/04 17:49:22  dicuccio
  302.  * Changed calling parameters for plugins - pass CPluginMessage instead of paired
  303.  * CPluginCommand/CPluginReply
  304.  *
  305.  * Revision 1.6  2003/11/03 17:41:19  dicuccio
  306.  * Fixed to match changes in BLAST API
  307.  *
  308.  * Revision 1.5  2003/10/23 16:21:28  dicuccio
  309.  * Remvoed dead code
  310.  *
  311.  * Revision 1.4  2003/10/16 21:52:07  jcherry
  312.  * Added filtering of neighbors
  313.  *
  314.  * ===========================================================================
  315.  */