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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: find_orfs.cpp,v $
  4.  * PRODUCTION Revision 1000.5  2004/06/01 20:54:59  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.26
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: find_orfs.cpp,v 1000.5 2004/06/01 20:54:59 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:  simple gbench plugin for finding ORFs
  37.  *
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include "find_orfs.hpp"
  41. #include <algo/sequence/make_cdr_prods.hpp>
  42. #include <algo/sequence/orf.hpp>
  43. #include <gui/core/plugin_utils.hpp>
  44. #include <gui/core/version.hpp>
  45. #include <gui/dialogs/col/multi_col_dlg.hpp>
  46. #include <gui/plugin/PluginCommandSet.hpp>
  47. #include <gui/plugin/PluginInfo.hpp>
  48. #include <gui/plugin/PluginReply.hpp>
  49. #include <gui/plugin/PluginRequest.hpp>
  50. #include <gui/plugin/PluginValueConstraint.hpp>
  51. #include <gui/objutils/utils.hpp>
  52. #include <objects/seqfeat/Genetic_code.hpp>
  53. #include <objects/seqfeat/Genetic_code_table.hpp>
  54. #include <objmgr/util/sequence.hpp>
  55. BEGIN_NCBI_SCOPE
  56. USING_SCOPE(objects);
  57. CAlgoPlugin_FindOrfs::~CAlgoPlugin_FindOrfs()
  58. {
  59. }
  60. // standard plugin announce bopilerplate
  61. void CAlgoPlugin_FindOrfs::GetInfo(CPluginInfo& info)
  62. {
  63.     info.Reset();
  64.     
  65.     // version info macro
  66.     info.SetInfo(CPluginVersion::eMajor, CPluginVersion::eMinor, 0,
  67.                  string(__DATE__) + " " + string(__TIME__),
  68.                  "CAlgoPlugin_FindOrfs", "Search/Find Open Reading Frames",
  69.                  "Find open reading frames in a DNA sequence",
  70.                  "");
  71.     // command info
  72.     CPluginCommandSet& cmds = info.SetCommands();
  73.     CPluginCommand&    args = cmds.AddAlgoCommand(eAlgoCommand_run);
  74.     args.AddArgument("locs", "Locations to evaluate",
  75.                      CSeq_loc::GetTypeInfo(),
  76.                      CPluginArg::TData::e_Array);
  77.     args.SetConstraint("locs",
  78.                        (*CPluginValueConstraint::CreateSeqMol(),
  79.                         CSeq_inst::eMol_na,
  80.                         CSeq_inst::eMol_dna,
  81.                         CSeq_inst::eMol_rna));
  82.     args.AddDefaultArgument("min_length_codons",
  83.                             "Minimum number of sense codons",
  84.                             CPluginArg::eInteger, "100");
  85.     // genetic code argument
  86.     const CGenetic_code_table& code_table = CGen_code_table::GetCodeTable();
  87.     const CGenetic_code_table::Tdata& codes = code_table.Get();
  88.     args.AddDefaultArgument("genetic_code", "Genetic code",
  89.                             CPluginArg::eString, codes.front()->GetName());
  90.     CPluginValueConstraint *code_list = CPluginValueConstraint::CreateSet();
  91.     ITERATE (CGenetic_code_table::Tdata, code, codes) {
  92.         code_list->SetSet().push_back((*code)->GetName());
  93.     }
  94.     args.SetConstraint("genetic_code", *code_list);
  95.                                        
  96. }
  97. void CAlgoPlugin_FindOrfs::RunCommand(CPluginMessage& msg)
  98. {
  99.     const CPluginCommand& args = msg.GetRequest().GetCommand();
  100.     CPluginReply& reply = msg.SetReply();
  101.     _TRACE("CAlgoPlugin_FindOrfs::RunCommand()");
  102.     
  103.     if ( !m_Dialog.get() ) {
  104.         m_Dialog.reset(new CMultiColDlg());
  105.         m_Dialog->SetWindowSize(500, 450);
  106.         m_Dialog->SetTitle("Open Reading Frames");
  107.         
  108.         m_Dialog->SetColumn(0, "Sequence", FL_ALIGN_LEFT, 0.5f);
  109.         m_Dialog->SetColumn(1, "Location", FL_ALIGN_LEFT, 0.5f);
  110.         m_Dialog->SetColumn(2, "Strand", FL_ALIGN_CENTER, 0.25f);
  111.         m_Dialog->SetColumn(3, "From", FL_ALIGN_CENTER, 0.5f);
  112.         m_Dialog->SetColumn(4, "To", FL_ALIGN_CENTER, 0.5f);
  113.         m_Dialog->SetColumn(5, "Sense Codons", FL_ALIGN_CENTER, 0.5f);
  114.     }
  115.     // clear any previous contents
  116.     m_Dialog->SetRows(0);
  117.     int row = 0;
  118.     plugin_args::TLocList locs;
  119.     GetArgValue(args["locs"], locs);
  120.     int min_length_codons = args["min_length_codons"].AsInteger();
  121.     string genetic_code_name = args["genetic_code"].AsString();
  122.     ITERATE (plugin_args::TLocList, iter, locs) {
  123.         const CSeq_loc&  loc = *iter->second;
  124.         const IDocument& doc = *iter->first;
  125.         // find the best ID for this bioseq
  126.         try {
  127.             CBioseq_Handle handle = doc.GetScope().GetBioseqHandle(loc);
  128.             // get sequence vector
  129.             CSeqVector vec =
  130.                 handle.GetSequenceView(loc,
  131.                                        CBioseq_Handle::eViewConstructed,
  132.                                        CBioseq_Handle::eCoding_Ncbi);
  133.             string& id_str  = m_Dialog->SetCell(row, 0);
  134.             string& loc_str = m_Dialog->SetCell(row, 1);
  135.             const CSeq_id& best_id =
  136.                 sequence::GetId(handle, sequence::eGetId_Best);
  137.             id_str.erase();
  138.             best_id.GetLabel(&id_str);
  139.             loc_str = CPluginUtils::GetLabel(loc, &doc.GetScope());
  140.             // place to store orfs
  141.             vector< CRef<CSeq_loc> > orfs;
  142.             // find some ORFs
  143.             COrf::FindOrfs(vec, orfs,
  144.                            min_length_codons * 3,
  145.                            x_DecodeGeneticCode(genetic_code_name));
  146.             // translate our locs to our parent location
  147.             NON_CONST_ITERATE (vector< CRef<CSeq_loc> >, iter, orfs) {
  148.                 (**iter).SetId(sequence::GetId(loc));
  149.                 *iter = CSeqUtils::RemapChildToParent(loc, **iter);
  150.             }
  151.             // make an annot
  152.             CRef<CSeq_id> this_id
  153.                 (const_cast<CSeq_id*>(&sequence::GetId(loc)));
  154.             CRef<CSeq_annot> annot =
  155.                 COrf::MakeCDSAnnot(orfs,
  156.                                    x_DecodeGeneticCode(genetic_code_name));
  157.             // add description to annot
  158.             annot->AddName("Open reading frames");
  159.             string comment =
  160.                 string("Open reading frames containing at least ") +
  161.                 NStr::IntToString(min_length_codons) +
  162.                 " sense codons using " + genetic_code_name +
  163.                 " genetic code";
  164.             annot->AddComment(comment);
  165.             // make protein sequences
  166.             CRef<CBioseq_set> product_set =
  167.                 CMakeCdrProds::MakeCdrProds(annot, handle);
  168.             reply.AddObject(doc, *product_set);
  169.             reply.AddObject(doc, *annot);
  170.             /**
  171.             CRef<CSeq_entry> new_entry(new CSeq_entry);
  172.             new_entry->SetSet(*product_set);
  173.             doc.GetScope().AddTopLevelSeqEntry(*new_entry);
  174.             **/
  175.             // attach annot to doc
  176.             //const_cast<IDocument&>(doc).AttachAnnot(*annot);
  177.             // in order to build dialog efficiently,
  178.             // pre-allocate one line for each ORF
  179.             m_Dialog->SetRows(row + orfs.size());
  180.             ITERATE (vector< CRef<CSeq_loc> >, loc_iter, orfs) {
  181.                 const CSeq_loc& orf = **loc_iter;
  182.                 //
  183.                 // add ORFs to dialog
  184.                 //
  185.                 ENa_strand strand = sequence::GetStrand(orf);
  186.                 if (strand == eNa_strand_minus) {
  187.                     m_Dialog->SetCell(row, 2) = "-";
  188.                 } else {
  189.                     m_Dialog->SetCell(row, 2) = "+";
  190.                 }
  191.                 m_Dialog->SetCell(row, 3)
  192.                     = NStr::IntToString(orf.GetTotalRange().GetFrom() + 1);
  193.                 m_Dialog->SetCell(row, 4)
  194.                     = NStr::IntToString(orf.GetTotalRange().GetTo() + 1);
  195.                 // ORF may or may not include a stop codon.
  196.                 // If it does, this must be subtracted
  197.                 // in computing the number of sense codons.
  198.                 int sense_codon_count = sequence::GetLength(orf);
  199.                 sense_codon_count /= 3;
  200.                 sense_codon_count -= 1;
  201.                 if ((strand == eNa_strand_plus   &&  orf.IsPartialRight())  ||
  202.                     (strand == eNa_strand_minus  &&  orf.IsPartialLeft())) {
  203.                     ++sense_codon_count;
  204.                 }
  205.                 m_Dialog->SetCell(row, 5)
  206.                     = NStr::IntToString(sense_codon_count);
  207.                 ++row;
  208.             }
  209.         }
  210.         catch (CException& e) {
  211.             LOG_POST(Error << "error processing location in ORF finder: "
  212.                      << e.what());
  213.             string str = CPluginUtils::GetLabel(loc, &doc.GetScope());
  214.             LOG_POST(Error << "Error processing location " << str);
  215.         }
  216.         catch (exception& e) {
  217.             LOG_POST(Error << "error processing location in ORF finder: "
  218.                      << e.what());
  219.             string str = CPluginUtils::GetLabel(loc, &doc.GetScope());
  220.             LOG_POST(Error << "Error processing location " << str);
  221.         }
  222. #ifndef _DEBUG
  223.         catch (...) {
  224.             string str = CPluginUtils::GetLabel(loc, &doc.GetScope());
  225.             LOG_POST(Error << "Error processing location " << str);
  226.         }
  227. #endif
  228.     }
  229.     // update all views
  230.     //CDocManager::UpdateAllViews();
  231.     //
  232.     // prepare our dialog box
  233.     //
  234.     m_Dialog->SetLabel(string("ORFs ") + NStr::IntToString(min_length_codons)
  235.                        + " codons or longer" 
  236.                        + " using " + genetic_code_name + " genetic code");
  237.     m_Dialog->Show();
  238.     reply.AddAction(CPluginReplyAction::e_Add_to_document);
  239.     reply.SetStatus(eMessageStatus_success);
  240. }
  241. // figure out the id of the genetic code the user wants
  242. int CAlgoPlugin_FindOrfs::x_DecodeGeneticCode(const string& s)
  243. {
  244.     const CGenetic_code_table& code_table = CGen_code_table::GetCodeTable();
  245.     const CGenetic_code_table::Tdata& codes = code_table.Get();
  246.     ITERATE (CGenetic_code_table::Tdata, code, codes) {
  247.         if ((*code)->GetName() == s) {
  248.             return (*code)->GetId();
  249.         }
  250.     }
  251.     // if we got here, nothing matched
  252.     NCBI_THROW(CException, eUnknown,
  253.                "CAlgoPlugin_FindOrfs: no genetic code matched " + s);
  254. }
  255. END_NCBI_SCOPE
  256. /*
  257.  * ===========================================================================
  258.  * $Log: find_orfs.cpp,v $
  259.  * Revision 1000.5  2004/06/01 20:54:59  gouriano
  260.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.26
  261.  *
  262.  * Revision 1.26  2004/05/21 22:27:46  gorelenk
  263.  * Added PCH ncbi_pch.hpp
  264.  *
  265.  * Revision 1.25  2004/05/03 13:05:42  dicuccio
  266.  * gui/utils --> gui/objutils where needed
  267.  *
  268.  * Revision 1.24  2004/03/05 17:35:07  dicuccio
  269.  * Use CGenetic_code_table typedefs to ease syntax.  Use sequence::GetId() instead
  270.  * of CSeq_id::GetStringDescr()
  271.  *
  272.  * Revision 1.23  2004/01/27 18:37:41  dicuccio
  273.  * Code clean-up.  Use standard names for plugins.  Removed unnecessary #includes
  274.  *
  275.  * Revision 1.22  2004/01/07 15:50:36  dicuccio
  276.  * Adjusted for API change in CPluginUtils::GetLabel().  Standardized exception
  277.  * reporting in algorithms.
  278.  *
  279.  * Revision 1.21  2003/11/24 15:45:26  dicuccio
  280.  * Renamed CVersion to CPluginVersion
  281.  *
  282.  * Revision 1.20  2003/11/18 17:48:36  dicuccio
  283.  * Added standard processing of return values
  284.  *
  285.  * Revision 1.19  2003/11/10 16:51:06  jcherry
  286.  * Added generation of protein sequences for orfs
  287.  *
  288.  * Revision 1.18  2003/11/06 20:12:12  dicuccio
  289.  * Cleaned up handling of USING_SCOPE - removed from all headers
  290.  *
  291.  * Revision 1.17  2003/11/04 17:49:22  dicuccio
  292.  * Changed calling parameters for plugins - pass CPluginMessage instead of paired
  293.  * CPluginCommand/CPluginReply
  294.  *
  295.  * Revision 1.16  2003/10/27 17:46:48  dicuccio
  296.  * Removed dead #includes
  297.  *
  298.  * Revision 1.15  2003/10/15 21:51:11  jcherry
  299.  * Don't set ids with MakeCDSAnnot; it doesn't work, and it would be
  300.  * redundant anyway.
  301.  *
  302.  * Revision 1.14  2003/10/15 13:40:26  dicuccio
  303.  * Mkae sure to set the 'id' for the seq-locs before calling RemapChildToParent()
  304.  *
  305.  * Revision 1.13  2003/10/14 16:24:37  dicuccio
  306.  * Correctly remap new feature locations through the parent location to the master
  307.  * sequence
  308.  *
  309.  * Revision 1.12  2003/10/07 13:47:00  dicuccio
  310.  * Renamed CPluginURL* to CPluginValue*
  311.  *
  312.  * Revision 1.11  2003/09/30 13:40:49  dicuccio
  313.  * Minor code clean-up: use container typedefs from ASN.1 generated classes
  314.  *
  315.  * Revision 1.10  2003/09/25 17:21:35  jcherry
  316.  * Added name to annot
  317.  *
  318.  * Revision 1.9  2003/09/04 19:27:53  jcherry
  319.  * Made an ORF include the stop codon, and marked certain ORFs as
  320.  * partial.  Put ability to construct a feature table into COrf.
  321.  *
  322.  * Revision 1.8  2003/09/04 14:05:24  dicuccio
  323.  * Use IDocument instead of CDocument
  324.  *
  325.  * Revision 1.7  2003/09/03 14:46:53  rsmith
  326.  * change namespace name from args to plugin_args to avoid clashes with variable names.
  327.  *
  328.  * Revision 1.6  2003/08/21 12:03:07  dicuccio
  329.  * Make use of new typedef in plugin_utils.hpp for argument values.
  330.  *
  331.  * Revision 1.5  2003/08/19 20:47:52  jcherry
  332.  * Use SetSet().pushback() rather than comma operator for adding
  333.  * constraints in loop (less bizarre-looking)
  334.  *
  335.  * Revision 1.4  2003/08/19 18:36:59  jcherry
  336.  * Allowed user to specify genetic code
  337.  *
  338.  * Revision 1.3  2003/08/18 19:24:15  jcherry
  339.  * Moved orf and seq_match to algo/sequence
  340.  *
  341.  * Revision 1.2  2003/08/18 18:01:58  jcherry
  342.  * Changed COrf::FindOrfs to produce a vector of CRef<CSeq_loc>.
  343.  * Added version of FindOrfs that takes a CSeqVector.
  344.  *
  345.  * Revision 1.1  2003/08/14 17:59:22  jcherry
  346.  * Initial version
  347.  *
  348.  * ===========================================================================
  349.  */