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

生物技术

开发平台:

C/C++

  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: reader_snp.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 19:41:43  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.12
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  9. /*  $Id: reader_snp.cpp,v 1000.1 2004/06/01 19:41:43 gouriano Exp $
  10.  * ===========================================================================
  11.  *                            PUBLIC DOMAIN NOTICE
  12.  *               National Center for Biotechnology Information
  13.  *
  14.  *  This software/database is a "United States Government Work" under the
  15.  *  terms of the United States Copyright Act.  It was written as part of
  16.  *  the author's official duties as a United States Government employee and
  17.  *  thus cannot be copyrighted.  This software/database is freely available
  18.  *  to the public for use. The National Library of Medicine and the U.S.
  19.  *  Government have not placed any restriction on its use or reproduction.
  20.  *
  21.  *  Although all reasonable efforts have been taken to ensure the accuracy
  22.  *  and reliability of the software and data, the NLM and the U.S.
  23.  *  Government do not and cannot warrant the performance or results that
  24.  *  may be obtained by using this software or data. The NLM and the U.S.
  25.  *  Government disclaim all warranties, express or implied, including
  26.  *  warranties of performance, merchantability or fitness for any particular
  27.  *  purpose.
  28.  *
  29.  *  Please cite the author in any work or product based on this material.
  30.  *
  31.  * ===========================================================================
  32.  *
  33.  *  Author:  Anton Butanaev, Eugene Vasilchenko
  34.  *
  35.  *  File Description: Data reader from ID1
  36.  *
  37.  */
  38. #include <ncbi_pch.hpp>
  39. #include <objtools/data_loaders/genbank/reader_snp.hpp>
  40. #include <objtools/data_loaders/genbank/reader.hpp>
  41. #include <objects/general/Object_id.hpp>
  42. #include <objects/general/User_object.hpp>
  43. #include <objects/general/User_field.hpp>
  44. #include <objects/general/Dbtag.hpp>
  45. #include <objects/seqloc/Seq_id.hpp>
  46. #include <objects/seqloc/Seq_point.hpp>
  47. #include <objects/seqloc/Seq_loc.hpp>
  48. #include <objects/seqfeat/Seq_feat.hpp>
  49. #include <objects/seqfeat/SeqFeatData.hpp>
  50. #include <objects/seqfeat/Imp_feat.hpp>
  51. #include <objects/seqfeat/Gb_qual.hpp>
  52. #include <objects/seqset/Seq_entry.hpp>
  53. #include <objects/seqset/Bioseq_set.hpp>
  54. #include <objects/seq/Seq_annot.hpp>
  55. #include <objmgr/objmgr_exception.hpp>
  56. #include <serial/objectinfo.hpp>
  57. #include <serial/objectiter.hpp>
  58. #include <serial/objectio.hpp>
  59. #include <serial/serial.hpp>
  60. #include <serial/objistr.hpp>
  61. #include <serial/objistrasnb.hpp>
  62. #include <serial/objostrasnb.hpp>
  63. #include <util/reader_writer.hpp>
  64. #include <algorithm>
  65. #include <numeric>
  66. // for debugging
  67. #include <serial/objostrasn.hpp>
  68. BEGIN_NCBI_SCOPE
  69. BEGIN_SCOPE(objects)
  70. class CSNP_Ftable_hook : public CReadChoiceVariantHook
  71. {
  72. public:
  73.     CSNP_Ftable_hook(CSeq_annot_SNP_Info& annot_snp_info)
  74.         : m_Seq_annot_SNP_Info(annot_snp_info)
  75.         {
  76.         }
  77.     void ReadChoiceVariant(CObjectIStream& in,
  78.                            const CObjectInfoCV& variant);
  79. private:
  80.     CSeq_annot_SNP_Info&   m_Seq_annot_SNP_Info;
  81. };
  82. class CSNP_Seq_feat_hook : public CReadContainerElementHook
  83. {
  84. public:
  85.     CSNP_Seq_feat_hook(CSeq_annot_SNP_Info& annot_snp_info,
  86.                        CSeq_annot::TData::TFtable& ftable);
  87.     ~CSNP_Seq_feat_hook(void);
  88.     void ReadContainerElement(CObjectIStream& in,
  89.                               const CObjectInfo& ftable);
  90. private:
  91.     CSeq_annot_SNP_Info&        m_Seq_annot_SNP_Info;
  92.     CSeq_annot::TData::TFtable& m_Ftable;
  93.     CRef<CSeq_feat>             m_Feat;
  94.     size_t                      m_Count[SSNP_Info::eSNP_Type_last];
  95. };
  96. void CSNP_Ftable_hook::ReadChoiceVariant(CObjectIStream& in,
  97.                                          const CObjectInfoCV& variant)
  98. {
  99.     CObjectInfo data_info = variant.GetChoiceObject();
  100.     CObjectInfo ftable_info = *variant;
  101.     CSeq_annot::TData& data = *CType<CSeq_annot::TData>::Get(data_info);
  102.     _ASSERT(ftable_info.GetObjectPtr() == static_cast<TConstObjectPtr>(&data.GetFtable()));
  103.     CSNP_Seq_feat_hook hook(m_Seq_annot_SNP_Info, data.SetFtable());
  104.     ftable_info.ReadContainer(in, hook);
  105. }
  106. static int s_GetEnvInt(const char* env, int def_val)
  107. {
  108.     const char* val = ::getenv(env);
  109.     if ( val ) {
  110.         try {
  111.             return NStr::StringToInt(val);
  112.         }
  113.         catch (...) {
  114.         }
  115.     }
  116.     return def_val;
  117. }
  118. static bool s_SNP_stat = s_GetEnvInt("GENBANK_SNP_TABLE_STAT", 0) > 0;
  119. CSNP_Seq_feat_hook::CSNP_Seq_feat_hook(CSeq_annot_SNP_Info& annot_snp_info,
  120.                                        CSeq_annot::TData::TFtable& ftable)
  121.     : m_Seq_annot_SNP_Info(annot_snp_info),
  122.       m_Ftable(ftable)
  123. {
  124.     fill(m_Count, m_Count+SSNP_Info::eSNP_Type_last, 0);
  125. }
  126. static size_t s_TotalCount[SSNP_Info::eSNP_Type_last] = { 0 };
  127. CSNP_Seq_feat_hook::~CSNP_Seq_feat_hook(void)
  128. {
  129.     if ( s_SNP_stat ) {
  130.         size_t total =
  131.             accumulate(m_Count, m_Count+SSNP_Info::eSNP_Type_last, 0);
  132.         NcbiCout << "CSeq_annot_SNP_Info statistic:n";
  133.         for ( size_t i = 0; i < SSNP_Info::eSNP_Type_last; ++i ) {
  134.             NcbiCout <<
  135.                 setw(40) << SSNP_Info::s_SNP_Type_Label[i] << ": " <<
  136.                 setw(6) << m_Count[i] << "  " <<
  137.                 setw(3) << int(m_Count[i]*100.0/total+.5) << "%n";
  138.             s_TotalCount[i] += m_Count[i];
  139.         }
  140.         NcbiCout << NcbiEndl;
  141.         total =
  142.             accumulate(s_TotalCount, s_TotalCount+SSNP_Info::eSNP_Type_last,0);
  143.         NcbiCout << "cumulative CSeq_annot_SNP_Info statistic:n";
  144.         for ( size_t i = 0; i < SSNP_Info::eSNP_Type_last; ++i ) {
  145.             NcbiCout <<
  146.                 setw(40) << SSNP_Info::s_SNP_Type_Label[i] << ": " <<
  147.                 setw(6) << s_TotalCount[i] << "  " <<
  148.                 setw(3) << int(s_TotalCount[i]*100.0/total+.5) << "%n";
  149.         }
  150.         NcbiCout << NcbiEndl;
  151.     }
  152. }
  153. void CSNP_Seq_feat_hook::ReadContainerElement(CObjectIStream& in,
  154.                                               const CObjectInfo& /*ftable*/)
  155. {
  156.     if ( !m_Feat ) {
  157.         m_Feat.Reset(new CSeq_feat);
  158.     }
  159.     in.ReadObject(&*m_Feat, m_Feat->GetTypeInfo());
  160.     SSNP_Info snp_info;
  161.     SSNP_Info::ESNP_Type snp_type =
  162.         snp_info.ParseSeq_feat(*m_Feat, m_Seq_annot_SNP_Info);
  163.     ++m_Count[snp_type];
  164.     if ( snp_type == SSNP_Info::eSNP_Simple ) {
  165.         m_Seq_annot_SNP_Info.x_AddSNP(snp_info);
  166.     }
  167.     else {
  168. #ifdef _DEBUG
  169.         static int dump_feature = s_GetEnvInt("GENBANK_SNP_TABLE_DUMP", 0);
  170.         if ( dump_feature ) {
  171.             --dump_feature;
  172.             NcbiCerr <<
  173.                 "CSNP_Seq_feat_hook::ReadContainerElement: complex SNP: " <<
  174.                 SSNP_Info::s_SNP_Type_Label[snp_type] << ":n" << 
  175.                 MSerial_AsnText << *m_Feat;
  176.         }
  177. #endif
  178.         m_Ftable.push_back(m_Feat);
  179.         m_Feat.Reset();
  180.     }
  181. }
  182. static void write_unsigned(CNcbiOstream& stream, unsigned n)
  183. {
  184.     stream.write(reinterpret_cast<const char*>(&n), sizeof(n));
  185. }
  186. static unsigned read_unsigned(CNcbiIstream& stream)
  187. {
  188.     unsigned n;
  189.     stream.read(reinterpret_cast<char*>(&n), sizeof(n));
  190.     return n;
  191. }
  192. static void write_size(CNcbiOstream& stream, unsigned size)
  193. {
  194.     // use ASN.1 binary like format
  195.     while ( size >= (1<<7) ) {
  196.         stream.put(char(size | (1<<7)));
  197.         size >>= 7;
  198.     }
  199.     stream.put(char(size));
  200. }
  201. static unsigned read_size(CNcbiIstream& stream)
  202. {
  203.     unsigned size = 0;
  204.     int shift = 0;
  205.     char c = char(1<<7);
  206.     while ( c & (1<<7) ) {
  207.         c = stream.get();
  208.         size |= (c & ((1<<7)-1)) << shift;
  209.         shift += 7;
  210.     }
  211.     return size;
  212. }
  213. void CIndexedStrings::StoreTo(CNcbiOstream& stream) const
  214. {
  215.     write_size(stream, m_Strings.size());
  216.     ITERATE ( TStrings, it, m_Strings ) {
  217.         unsigned size = it->size();
  218.         write_size(stream, size);
  219.         stream.write(it->data(), size);
  220.     }
  221. }
  222. void CIndexedStrings::LoadFrom(CNcbiIstream& stream,
  223.                                size_t max_index,
  224.                                size_t max_length)
  225. {
  226.     Clear();
  227.     unsigned count = read_size(stream);
  228.     if ( !stream || (count > unsigned(max_index+1)) ) {
  229.         NCBI_THROW(CLoaderException, eLoaderFailed,
  230.                    "Bad format of SNP table");
  231.     }
  232.     m_Strings.resize(count);
  233.     AutoPtr<char, ArrayDeleter<char> > buf(new char[max_length]);
  234.     NON_CONST_ITERATE ( TStrings, it, m_Strings ) {
  235.         unsigned size = read_size(stream);
  236.         if ( !stream || (size > max_length) ) {
  237.             Clear();
  238.             NCBI_THROW(CLoaderException, eLoaderFailed,
  239.                        "Bad format of SNP table");
  240.         }
  241.         stream.read(buf.get(), size);
  242.         if ( !stream ) {
  243.             Clear();
  244.             NCBI_THROW(CLoaderException, eLoaderFailed,
  245.                        "Bad format of SNP table");
  246.         }
  247.         it->assign(buf.get(), buf.get()+size);
  248.     }
  249. }
  250. void CSeq_annot_SNP_Info_Reader::Parse(CObjectIStream& in,
  251.                                        CSeq_annot_SNP_Info& snp_info)
  252. {
  253.     snp_info.m_Seq_annot.Reset(new CSeq_annot); // Seq-annot object
  254.     CReader::SetSNPReadHooks(in);
  255.     if ( CReader::TrySNPTable() ) { // set SNP hook
  256.         CObjectTypeInfo type = CType<CSeq_annot::TData>();
  257.         CObjectTypeInfoVI ftable = type.FindVariant("ftable");
  258.         ftable.SetLocalReadHook(in, new CSNP_Ftable_hook(snp_info));
  259.     }
  260.     
  261.     in >> *snp_info.m_Seq_annot;
  262.     // we don't need index maps anymore
  263.     snp_info.m_Comments.ClearIndices();
  264.     snp_info.m_Alleles.ClearIndices();
  265.     sort(snp_info.m_SNP_Set.begin(), snp_info.m_SNP_Set.end());
  266. }
  267. static const unsigned MAGIC = 0x12340002;
  268. void CSeq_annot_SNP_Info_Reader::Write(CNcbiOstream& stream,
  269.                                        const CSeq_annot_SNP_Info& snp_info)
  270. {
  271.     // header
  272.     write_unsigned(stream, MAGIC);
  273.     write_unsigned(stream, snp_info.GetGi());
  274.     // strings
  275.     snp_info.m_Comments.StoreTo(stream);
  276.     snp_info.m_Alleles.StoreTo(stream);
  277.     // simple SNPs
  278.     unsigned count = snp_info.m_SNP_Set.size();
  279.     write_size(stream, count);
  280.     stream.write(reinterpret_cast<const char*>(&snp_info.m_SNP_Set[0]),
  281.                  count*sizeof(SSNP_Info));
  282.     // complex SNPs
  283.     CObjectOStreamAsnBinary obj_stream(stream);
  284.     obj_stream << *snp_info.m_Seq_annot;
  285. }
  286. void CSeq_annot_SNP_Info_Reader::Read(CNcbiIstream& stream,
  287.                                       CSeq_annot_SNP_Info& snp_info)
  288. {
  289.     snp_info.Reset();
  290.     // header
  291.     unsigned magic = read_unsigned(stream);
  292.     if ( !stream || magic != MAGIC ) {
  293.         NCBI_THROW(CLoaderException, eLoaderFailed,
  294.                    "Incompatible version of SNP table");
  295.     }
  296.     snp_info.x_SetGi(read_unsigned(stream));
  297.     // strings
  298.     snp_info.m_Comments.LoadFrom(stream,
  299.                                  SSNP_Info::kMax_CommentIndex,
  300.                                  SSNP_Info::kMax_CommentLength);
  301.     snp_info.m_Alleles.LoadFrom(stream,
  302.                                 SSNP_Info::kMax_AlleleIndex,
  303.                                 SSNP_Info::kMax_AlleleLength);
  304.     // simple SNPs
  305.     unsigned count = read_size(stream);
  306.     if ( stream ) {
  307.         snp_info.m_SNP_Set.resize(count);
  308.         stream.read(reinterpret_cast<char*>(&snp_info.m_SNP_Set[0]),
  309.                     count*sizeof(SSNP_Info));
  310.     }
  311.     size_t comments_size = snp_info.m_Comments.GetSize();
  312.     size_t alleles_size = snp_info.m_Alleles.GetSize();
  313.     ITERATE ( CSeq_annot_SNP_Info::TSNP_Set, it, snp_info.m_SNP_Set ) {
  314.         size_t index = it->m_CommentIndex;
  315.         if ( index != SSNP_Info::kNo_CommentIndex &&
  316.              index >= comments_size ) {
  317.             snp_info.Reset();
  318.             NCBI_THROW(CLoaderException, eLoaderFailed,
  319.                        "Bad format of SNP table");
  320.         }
  321.         for ( size_t i = 0; i < SSNP_Info::kMax_AllelesCount; ++i ) {
  322.             index = it->m_AllelesIndices[i];
  323.             if ( index != SSNP_Info::kNo_AlleleIndex &&
  324.                  index >= alleles_size ) {
  325.                 snp_info.Reset();
  326.                 NCBI_THROW(CLoaderException, eLoaderFailed,
  327.                            "Bad format of SNP table");
  328.             }
  329.         }
  330.     }
  331.     // complex SNPs
  332.     CObjectIStreamAsnBinary obj_stream(stream);
  333.     if ( !snp_info.m_Seq_annot ) {
  334.         snp_info.m_Seq_annot.Reset(new CSeq_annot);
  335.     }
  336.     obj_stream >> *snp_info.m_Seq_annot;
  337.     if ( !stream ) {
  338.         snp_info.m_Seq_annot.Reset();
  339.         NCBI_THROW(CLoaderException, eLoaderFailed,
  340.                    "Bad format of SNP table");
  341.     }
  342. }
  343. END_SCOPE(objects)
  344. END_NCBI_SCOPE
  345. /*
  346.  * $Log: reader_snp.cpp,v $
  347.  * Revision 1000.1  2004/06/01 19:41:43  gouriano
  348.  * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.12
  349.  *
  350.  * Revision 1.12  2004/05/21 21:42:52  gorelenk
  351.  * Added PCH ncbi_pch.hpp
  352.  *
  353.  * Revision 1.11  2004/03/16 16:04:20  vasilche
  354.  * Removed conversion warning
  355.  *
  356.  * Revision 1.10  2004/02/06 16:13:19  vasilche
  357.  * Added parsing "replace" as a synonym of "allele" in SNP qualifiers.
  358.  * More compact format of SNP table in cache. SNP table version increased.
  359.  * Fixed null pointer exception when SNP features are loaded from cache.
  360.  *
  361.  * Revision 1.9  2004/01/13 16:55:55  vasilche
  362.  * CReader, CSeqref and some more classes moved from xobjmgr to separate lib.
  363.  * Headers moved from include/objmgr to include/objtools/data_loaders/genbank.
  364.  *
  365.  * Revision 1.8  2003/10/23 13:47:56  vasilche
  366.  * Fixed Reset() method: m_Gi should be -1.
  367.  *
  368.  * Revision 1.7  2003/10/21 16:29:13  vasilche
  369.  * Added check for errors in SNP table loaded from cache.
  370.  *
  371.  * Revision 1.6  2003/10/21 15:21:21  vasilche
  372.  * Avoid use of non-constant array sizes of stack arrays.
  373.  *
  374.  * Revision 1.5  2003/10/21 14:27:35  vasilche
  375.  * Added caching of gi -> sat,satkey,version resolution.
  376.  * SNP blobs are stored in cache in preprocessed format (platform dependent).
  377.  * Limit number of connections to GenBank servers.
  378.  * Added collection of ID1 loader statistics.
  379.  *
  380.  * Revision 1.4  2003/09/30 16:22:03  vasilche
  381.  * Updated internal object manager classes to be able to load ID2 data.
  382.  * SNP blobs are loaded as ID2 split blobs - readers convert them automatically.
  383.  * Scope caches results of requests for data to data loaders.
  384.  * Optimized CSeq_id_Handle for gis.
  385.  * Optimized bioseq lookup in scope.
  386.  * Reduced object allocations in annotation iterators.
  387.  * CScope is allowed to be destroyed before other objects using this scope are
  388.  * deleted (feature iterators, bioseq handles etc).
  389.  * Optimized lookup for matching Seq-ids in CSeq_id_Mapper.
  390.  * Added 'adaptive' option to objmgr_demo application.
  391.  *
  392.  * Revision 1.3  2003/08/19 18:35:21  vasilche
  393.  * CPackString classes were moved to SERIAL library.
  394.  *
  395.  * Revision 1.2  2003/08/15 19:19:16  vasilche
  396.  * Fixed memory leak in string packing hooks.
  397.  * Fixed processing of 'partial' flag of features.
  398.  * Allow table packing of non-point SNP.
  399.  * Allow table packing of SNP with long alleles.
  400.  *
  401.  * Revision 1.1  2003/08/14 20:05:19  vasilche
  402.  * Simple SNP features are stored as table internally.
  403.  * They are recreated when needed using CFeat_CI.
  404.  *
  405.  */