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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: align_5prime.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 20:54:04  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: align_5prime.cpp,v 1000.1 2004/06/01 20:54:04 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:  Mike DiCuccio
  35.  *
  36.  * File Description:
  37.  *    gbench plugin for aligning to neighbors
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include "align_5prime.hpp"
  41. #include "blast_util.hpp"
  42. #include <algo/blast/api/bl2seq.hpp>
  43. #include <gui/core/doc_manager.hpp>
  44. #include <gui/core/plugin_utils.hpp>
  45. #include <gui/core/version.hpp>
  46. #include <gui/plugin/AlgoCommand.hpp>
  47. #include <gui/plugin/PluginCommand.hpp>
  48. #include <gui/plugin/PluginCommandSet.hpp>
  49. #include <gui/plugin/PluginInfo.hpp>
  50. #include <gui/plugin/PluginReply.hpp>
  51. #include <gui/plugin/PluginRequest.hpp>
  52. #include <gui/plugin/PluginValueConstraint.hpp>
  53. #include <objects/seqalign/Seq_align.hpp>
  54. #include <objects/seqfeat/Seq_feat.hpp>
  55. #include <objects/seqloc/Seq_loc.hpp>
  56. #include <objmgr/util/feature.hpp>
  57. #include <objmgr/util/sequence.hpp>
  58. BEGIN_NCBI_SCOPE
  59. USING_SCOPE(ncbi::objects);
  60. void CAlgoPlugin_Align5Prime::GetInfo(CPluginInfo& info)
  61. {
  62.     info.Reset();
  63.     // version info macro
  64.     info.SetInfo(CPluginVersion::eMajor, CPluginVersion::eMinor, 0,
  65.                  string(__DATE__) + " " + string(__TIME__),
  66.                  "CAlgoPlugin_Align5Prime",
  67.                  "Alignments/Align 5' region of genes",
  68.                  "Create an alignment of the 5' flanking regions of genes",
  69.                  "");
  70.     // command info
  71.     CPluginCommandSet& cmds = info.SetCommands();
  72.     CPluginCommand&    args = cmds.AddAlgoCommand(eAlgoCommand_run);
  73.     args.AddArgument("genes", "Genes to align",
  74.                      CSeq_feat::GetTypeInfo(),
  75.                      CPluginArg::TData::e_Array);
  76.     args.SetConstraint("genes",
  77.                        (*CPluginValueConstraint::CreateFeatSubtype(),
  78.                         CSeqFeatData::eSubtype_gene));
  79.     args.AddDefaultArgument("offset", "Alignment offset",
  80.                             CPluginArg::eInteger, "10000");
  81.     args.AddDefaultArgument("length", "Alignment length",
  82.                             CPluginArg::eInteger, "10000");
  83.     CBlastUtils::AddBlastArgs(args, blast::eBlastn);
  84. }
  85. void CAlgoPlugin_Align5Prime::RunCommand(CPluginMessage& msg)
  86. {
  87.     const CPluginCommand& args = msg.GetRequest().GetCommand();
  88.     CPluginReply&         reply = msg.SetReply();
  89.     reply.SetStatus(eMessageStatus_failed);
  90.     plugin_args::TFeatList feats;
  91.     GetArgValue(args["genes"], feats);
  92.     TSeqPos offs = args["offset"].AsInteger();
  93.     TSeqPos len  = args["offset"].AsInteger();
  94.     // scan our features.  We must make sure that we get the appropriate
  95.     // upstream location - some features may be too long
  96.     blast::TSeqLocVector targets;
  97.     ITERATE (plugin_args::TFeatList, iter, feats) {
  98.         const CSeq_feat& feat = *iter->second;
  99.         const IDocument& doc  = *iter->first;
  100.         CRef<CSeq_loc> loc(new CSeq_loc());
  101.         ENa_strand strand = sequence::GetStrand(feat.GetLocation());
  102.         TSeqRange range = feat.GetLocation().GetTotalRange();
  103.         TSeqPos from;
  104.         TSeqPos to;
  105.         if (strand == eNa_strand_minus) {
  106.             from = range.GetFrom() + offs;
  107.             // make sure we're bounded by the sequence upstream
  108.             CBioseq_Handle handle =
  109.                 doc.GetScope().GetBioseqHandle(feat.GetLocation());
  110.             if (from > handle.GetBioseqLength()) {
  111.                 from = handle.GetBioseqLength();
  112.             }
  113.             to = from - len;
  114.             if (len > from) {
  115.                 to = 0;
  116.             }
  117.         } else {
  118.             // assume '+' strand
  119.             from = range.GetFrom() - offs;
  120.             // make sure we're bounded by the sequence upstream
  121.             if (range.GetFrom() < offs) {
  122.                 from = 0;
  123.             }
  124.             to = from + len;
  125.             CBioseq_Handle handle =
  126.                 doc.GetScope().GetBioseqHandle(feat.GetLocation());
  127.             if (to > handle.GetBioseqLength()) {
  128.                 to = handle.GetBioseqLength();
  129.             }
  130.         }
  131.         loc->SetInt().SetFrom(from);
  132.         loc->SetInt().SetTo  (to);
  133.         loc->SetInt().SetStrand(strand);
  134.         loc->SetId(sequence::GetId(feat.GetLocation()));
  135.         blast::SSeqLoc blast_loc(loc, &doc.GetScope());
  136.         targets.push_back(blast_loc);
  137.     }
  138.     // find a long one to be the query
  139.     blast::TSeqLocVector::iterator iter = targets.begin();
  140.     blast::TSeqLocVector::iterator best = targets.begin();
  141.     ++iter;
  142.     for ( ;  iter != targets.end();  ++iter) {
  143.         if (iter->seqloc->GetTotalRange().GetLength() >
  144.             best->seqloc->GetTotalRange().GetLength()) {
  145.             best = iter;
  146.         }
  147.     }
  148.     blast::SSeqLoc query(*best);
  149.     targets.erase(best);
  150.     blast::CBl2Seq blaster(query, targets, blast::eBlastn);
  151.     CBlastUtils::ArgsToBlastOptions(args, blaster.SetOptions());
  152.     // run blast
  153. blast::TSeqAlignVector aligns = blaster.Run();
  154.     // make an annotation
  155.     CRef<CSeq_annot> annot(new CSeq_annot());
  156. ITERATE (blast::TSeqAlignVector, iter, aligns) {
  157. annot->SetData().SetAlign()
  158. .insert(annot->SetData().SetAlign().end(),
  159. (*iter)->Get().begin(), (*iter)->Get().end());
  160. }
  161.     annot->AddName("Upstream flanking alignment");
  162.     //
  163.     // pass back to the system.  We may use the same scope and just attach,
  164.     // if that is appropriate
  165.     //
  166.     CConstRef<IDocument> doc_ref(feats.begin()->first);
  167.     ITERATE (plugin_args::TFeatList, iter, feats) {
  168.         if ( !doc_ref ) {
  169.             doc_ref.Reset(iter->first);
  170.         } else if (iter->first != doc_ref) {
  171.             doc_ref.Reset();
  172.             break;
  173.         }
  174.     }
  175.     if ( !doc_ref ) {
  176.         //
  177.         // features come from different documents
  178.         // create a new one to handle the results
  179.         //
  180.         CRef<CScope> new_scope(new CScope(CDocManager::GetObjectManager()));
  181.         ITERATE (plugin_args::TFeatList, iter, feats) {
  182.             new_scope->AddScope(iter->first->GetScope());
  183.         }
  184.         doc_ref.Reset(CDocManager::CreateDocument(*new_scope, *annot));
  185.     } else {
  186.         reply.AddAction(CPluginReplyAction::e_Add_to_document);
  187.     }
  188.     // standard reply handling
  189.     reply.AddObject(*doc_ref, *annot);
  190.     reply.SetStatus(eMessageStatus_success);
  191.     // go ahead and create an alignment view
  192.     CRef<CSelectionBuffer> buf(new CSelectionBuffer());
  193.     buf->AddSelection(doc_ref, *annot);
  194.     CPluginUtils::CallPlugin("CAlnMultiView",
  195.                              eViewCommand_new_view, *buf);
  196. }
  197. END_NCBI_SCOPE
  198. /*
  199.  * ===========================================================================
  200.  * $Log: align_5prime.cpp,v $
  201.  * Revision 1000.1  2004/06/01 20:54:04  gouriano
  202.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.2
  203.  *
  204.  * Revision 1.2  2004/05/21 22:27:46  gorelenk
  205.  * Added PCH ncbi_pch.hpp
  206.  *
  207.  * Revision 1.1  2004/01/27 18:39:01  dicuccio
  208.  * Initial revision
  209.  *
  210.  * Revision 1.11  2004/01/13 20:36:29  dicuccio
  211.  * ===========================================================================
  212.  */