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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: alignment_smear.cpp,v $
  4.  * PRODUCTION Revision 1000.0  2004/06/01 21:19:53  gouriano
  5.  * PRODUCTION PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.2
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: alignment_smear.cpp,v 1000.0 2004/06/01 21:19:53 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: Robert Smith
  35.  *
  36.  * File Description:
  37.  *
  38.  */
  39. #include <ncbi_pch.hpp>
  40. #include <gui/objutils/alignment_smear.hpp>
  41. #include <objtools/alnmgr/alnmap.hpp>
  42. #include <objmgr/bioseq_handle.hpp>
  43. #include <objmgr/annot_selector.hpp>
  44. #include <objects/seq/Seq_annot.hpp>
  45. #include <objects/seq/Annot_descr.hpp>
  46. #include <objects/seq/Annotdesc.hpp>
  47. #include <objects/general/User_object.hpp>
  48. #include <objects/general/User_field.hpp>
  49. #include <objects/general/Object_id.hpp>
  50. BEGIN_NCBI_SCOPE
  51. USING_SCOPE(objects);
  52. CAlignmentSmear::CAlignmentSmear(
  53.         const objects::CBioseq_Handle& handle, 
  54.         TSeqPos start, 
  55.         TSeqPos stop, 
  56.         TSeqPos window,
  57.         EAlignSmearStrand strand_type
  58.     )
  59.     : m_BioseqHandle(handle),
  60.     m_AccumSeg(start, stop, window),
  61.     m_AccumGap(start, stop, window),
  62.     m_StrandType(strand_type)
  63. {
  64. }
  65. /// smear all the alignments in this annotation.
  66. void CAlignmentSmear::AddAnnot(const CSeq_annot& seq_annot)
  67. {
  68.     objects::SAnnotSelector sel;
  69.     sel.SetLimitSeqAnnot(&seq_annot)
  70.        .SetSortOrder(SAnnotSelector::eSortOrder_None);
  71.     AddAlignments(sel);
  72.     
  73.     // now if there is a chance we would put more than one annotation
  74.     // in a smear we need to join their names together here, or something else?
  75.     SetLabel(x_GetAnnotName(seq_annot));
  76. }
  77. #if 1
  78. /// Smear all the alignments matchec by this selector on my bioseq.
  79. void CAlignmentSmear::AddAlignments(const objects::SAnnotSelector& sel)
  80. {
  81.     TSeqPos start(m_AccumSeg.GetStart());
  82.     TSeqPos stop(m_AccumSeg.GetStop());
  83.       
  84.     
  85.    // grab a alignment iterator for our range
  86.     CAlign_CI align_iter(m_BioseqHandle, start, stop, sel);
  87.     for (int ai = 0; align_iter; ++align_iter, ++ai) {
  88.         const CSeq_align& align = *align_iter;
  89.         
  90.         // convert to an AlnMap
  91.         CRef< CAlnMap> aln_map;
  92.         if (align.GetSegs().IsStd()) {
  93.             continue;   // TODO: Make  StdSeqs work.
  94.             CRef<CSeq_align> ds_align = align.CreateDensegFromStdseg();
  95.             aln_map =  new CAlnMap( ds_align->GetSegs().GetDenseg());
  96.         } else if (align.GetSegs().IsDenseg()) {
  97.             aln_map = new CAlnMap(align.GetSegs().GetDenseg());
  98.         }
  99.         AddAlignment(*aln_map);
  100.     }
  101.     MaskGaps();
  102. }
  103. #else
  104. /* same thing but we group alignments by their id */
  105. /// Smear all the alignments matchec by this selector on my bioseq.
  106. void CAlignmentSmear::AddAlignments(const objects::SAnnotSelector& sel)
  107. {
  108.     TSeqPos start(m_AccumSeg.GetStart());
  109.     TSeqPos stop(m_AccumSeg.GetStop());
  110.       
  111.     typedef map<string, CRef<CAlnMix> > TMixes;
  112.     TMixes  mixes;
  113.     
  114.    // grab a alignment iterator for our range
  115.     CAlign_CI align_iter(m_BioseqHandle, start, stop, sel);
  116.     for (int ai = 0; align_iter; ++align_iter, ++ai) {
  117.         const CSeq_align& align = *align_iter;
  118.         // convert to an AlnMap
  119.         CRef< CAlnMap> aln_map;
  120.         if (align.GetSegs().IsStd()) {
  121.             CRef<CSeq_align> ds_align = align.CreateDensegFromStdseg();
  122.             aln_map =  new CAlnMap( ds_align->GetSegs().GetDenseg());
  123.         } else if (align.GetSegs().IsDenseg()) {
  124.             aln_map = new CAlnMap(align.GetSegs().GetDenseg());
  125.         }
  126.         // which row of the alignment is the other sequence on?
  127.         CAlnMap::TNumrow row = m_BioseqHandle.IsSynonym( 
  128.                             aln_map->GetSeqId(0)) ? 1 : 0;
  129.                             
  130.         // What is that sequences label?
  131.         string  align_label;
  132.         aln_map->GetSeqId(row).GetLabel(&align_label);        
  133.         // Mix the alignments based on their label.
  134.         TMixes::const_iterator mit = mixes.find(align_label);
  135.         if (mit == mixes.end()) {
  136.             CRef<CAlnMix> mix_ref(new CAlnMix(m_BioseqHandle.GetScope()));
  137.             mixes[align_label] = mix_ref;
  138.         } 
  139.         mixes[align_label]->Add(aln_map->GetDenseg());
  140.     }
  141.     
  142.     // Now smear the mixed alignments.
  143.     NON_CONST_ITERATE(TMixes, mit, mixes)
  144.     {
  145.         CAlnMix & aln_mix(*(mit->second));
  146.         aln_mix.Merge(CAlnMix::fGen2EST |
  147.                           CAlnMix::fTryOtherMethodOnFail |
  148.                           CAlnMix::fGapJoin);
  149.          
  150.         // convert to a CAlnMap
  151.         CAlnMap aln_map(aln_mix.GetDenseg());
  152.         AddAlignment(aln_map);
  153.     }
  154.     
  155.     // take out gaps in the middle of blocks.
  156.     MaskGaps();
  157. }
  158. #endif
  159. void CAlignmentSmear::AddAlignment(CAlnMap& aln_map)
  160. {
  161.     // check number of dimensions
  162.     if ( aln_map.GetNumRows() != 2) {
  163.         // cout << "alignment does not have 2 rows." << endl;
  164.         return;
  165.     }
  166.     // find out what row our bioseq is on.        
  167.     CAlnMap::TNumrow bs_row = m_BioseqHandle.IsSynonym( 
  168.                             aln_map.GetSeqId(0)) ? 0 : 1;
  169.     _ASSERT( bs_row == 0  ||  bs_row == 1 );
  170.     CAlnMap::TNumrow other_row = 1 - bs_row; // works since bs_row is always 0 or 1.
  171.     // Are we only doing one strand?
  172.     if (m_StrandType == eSmearStrand_Pos  && aln_map.IsNegativeStrand(bs_row) )
  173.         return;
  174.     if (m_StrandType == eSmearStrand_Neg  && aln_map.IsPositiveStrand(bs_row) )
  175.         return;
  176.     
  177.     CAlnMap::TSignedRange range(aln_map.GetAlnStart(),
  178.                                 aln_map.GetAlnStop());
  179.     CRef<CAlnMap::CAlnChunkVec> aln_chunks
  180.             (aln_map.GetAlnChunks(other_row, range,
  181.                        CAlnMap::fSeqOnly )); 
  182.                        // | CAlnMap::fChunkSameAsSeg));
  183.                        
  184.     // Smear all the segments.
  185.     for (size_t i = 0;  i < aln_chunks->size();  ++i) {
  186.         CConstRef<CAlnMap::CAlnChunk> chunk((*aln_chunks)[i]);
  187.         
  188.         TSeqPos start = chunk->GetRange().GetFrom();
  189.         start = aln_map.GetSeqPosFromSeqPos(bs_row, other_row, start);
  190.         TSeqPos stop = chunk->GetRange().GetTo();
  191.         stop = aln_map.GetSeqPosFromSeqPos(bs_row, other_row, stop);
  192.         
  193.         m_AccumSeg.AddRange(TSeqRange(start, stop), 1);
  194.     }
  195.     
  196.     // Smear the gaps.
  197.     for (size_t i = 1;  i < aln_chunks->size();  ++i) {
  198.         CConstRef<CAlnMap::CAlnChunk> prev_chunk((*aln_chunks)[i-1]);
  199.         CConstRef<CAlnMap::CAlnChunk> chunk((*aln_chunks)[i]);
  200.         TSeqPos start = prev_chunk->GetRange().GetTo();
  201.         start = aln_map.GetSeqPosFromSeqPos(bs_row, other_row, start);
  202.         TSeqPos stop = chunk->GetRange().GetFrom();
  203.         stop = aln_map.GetSeqPosFromSeqPos(bs_row, other_row, stop);
  204.         
  205.         m_AccumGap.AddRange(TSeqRange(start, stop));
  206.     }
  207. }
  208. // erase gaps where there are segments.
  209. void CAlignmentSmear::MaskGaps()
  210. {
  211.     TSegMap::iterator seg_it = m_AccumSeg.begin();
  212.     TGapMap::iterator gap_it = m_AccumGap.begin();
  213.     for (; gap_it != m_AccumGap.end(); ++gap_it, ++seg_it  ) {
  214.         if (*gap_it > 0   &&  *seg_it > 0) {
  215.             *gap_it = 0;
  216.         }
  217.     }
  218.     
  219. }
  220. string CAlignmentSmear::GetLabel() const
  221. {
  222.     return m_Label;
  223. }
  224. void CAlignmentSmear::SetLabel(const string& label)
  225. {
  226.     m_Label = label;
  227. }
  228. bool CAlignmentSmear::SeparateStrands(const objects::CSeq_annot& seq_annot)
  229. {
  230.     string name = x_GetAnnotName(seq_annot);
  231.     
  232.     if (name == "BLASTX - swissprot"  || 
  233.         name == "BLASTN - mrna" ) {
  234.         return true;   
  235.     }
  236.     return false;
  237. }
  238. static const string s_BlastFieldKey("Blast Type");
  239. static const string s_HistFieldKey("Hist Seqalign");
  240. string CAlignmentSmear::x_GetAnnotName(const objects::CSeq_annot& seq_annot)
  241. {
  242.     string label;
  243.     string annot_name;
  244.     string blast_type;
  245.     string hist_type;
  246.     
  247.     if (seq_annot.IsSetDesc()) {
  248.         ITERATE ( CSeq_annot::TDesc::Tdata, it,
  249.                   seq_annot.GetDesc().Get() ) {
  250.             const CAnnotdesc& desc = **it;
  251.             if (desc.Which() == CAnnotdesc::e_Name ) {
  252.                 annot_name = desc.GetName();
  253.             }
  254.             if (desc.Which() == CAnnotdesc::e_User) {
  255.                 const CUser_object& user_obj = desc.GetUser();
  256.             
  257.                 if (user_obj.CanGetType()  &&
  258.                     user_obj.GetType().Which() == CObject_id::e_Str ) {
  259.                 
  260.                     if (user_obj.GetType().GetStr()  ==  s_BlastFieldKey) {
  261.                         const CUser_field& user_field = * user_obj.GetData().front();
  262.                         
  263.                         if (user_field.CanGetData()  &&  
  264.                             user_field.GetData().Which() == CUser_field::TData::e_Int &&
  265.                             user_field.CanGetLabel()) {
  266.                             const CObject_id& id = user_field.GetLabel();
  267.                             
  268.                             if (id.Which() == CObject_id::e_Str ) {
  269.                                 blast_type = id.GetStr();
  270.                             }
  271.                         }
  272.                     }
  273.                 
  274.                     if (user_obj.GetType().GetStr()  ==  s_HistFieldKey  ) {
  275.                         const CUser_field& user_field = * user_obj.GetData().front();
  276.                         if (user_field.CanGetData()  &&  
  277.                             user_field.GetData().Which() == CUser_field::TData::e_Bool &&
  278.                             user_field.GetData().GetBool()  &&
  279.                             user_field.CanGetLabel()) {
  280.                             const CObject_id& id = user_field.GetLabel();
  281.                             if (id.Which() == CObject_id::e_Str ) {
  282.                                 hist_type = id.GetStr();
  283.                             }
  284.                         }
  285.                     }
  286.                 }
  287.             }
  288.         }
  289.         if ( ! annot_name.empty()) {
  290.             label = annot_name;
  291.         } else if ( ! blast_type.empty() ) {
  292.             label = blast_type;
  293.         } else if ( ! hist_type.empty() ) {
  294.             label = hist_type;
  295.         }
  296.     } 
  297.     return label;
  298. }
  299. END_NCBI_SCOPE
  300. /*
  301.  * ===========================================================================
  302.  * $Log: alignment_smear.cpp,v $
  303.  * Revision 1000.0  2004/06/01 21:19:53  gouriano
  304.  * PRODUCTION: IMPORTED [GCC34_MSVC7] Dev-tree R1.2
  305.  *
  306.  * Revision 1.2  2004/05/21 22:27:43  gorelenk
  307.  * Added PCH ncbi_pch.hpp
  308.  *
  309.  * Revision 1.1  2004/04/30 11:48:15  dicuccio
  310.  * Initial commit - split out from src/gui/utils
  311.  *
  312.  * Revision 1.2  2004/04/07 15:22:39  rsmith
  313.  * new (not yet used) implementation of AddAlignments
  314.  *
  315.  * Revision 1.1  2004/03/22 16:40:55  rsmith
  316.  * initial checkin
  317.  *
  318.  *
  319.  * ===========================================================================
  320.  */