开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
gff_reader.cpp
资源名称：ncbi_cxx_win--Jul_08_2004.zip　[点击查看]
上传用户：yhdzpy8989
上传日期：2007-06-13
资源大小：13604k
文件大小：24k
源码类别：

生物技术

开发平台：

C/C++

gff_reader.cpp：源码内容
  1. /*
  2.  * ===========================================================================
  3.  * PRODUCTION $Log: gff_reader.cpp,v $
  4.  * PRODUCTION Revision 1000.1  2004/06/01 19:46:20  gouriano
  5.  * PRODUCTION PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  6.  * PRODUCTION
  7.  * ===========================================================================
  8.  */
  10. /*  $Id: gff_reader.cpp,v 1000.1 2004/06/01 19:46:20 gouriano Exp $
  11. * ===========================================================================
  12. *
  13. *                            PUBLIC DOMAIN NOTICE
  14. *               National Center for Biotechnology Information
  15. *
  16. *  This software/database is a "United States Government Work" under the
  17. *  terms of the United States Copyright Act.  It was written as part of
  18. *  the author's official duties as a United States Government employee and
  19. *  thus cannot be copyrighted.  This software/database is freely available
  20. *  to the public for use. The National Library of Medicine and the U.S.
  21. *  Government have not placed any restriction on its use or reproduction.
  22. *
  23. *  Although all reasonable efforts have been taken to ensure the accuracy
  24. *  and reliability of the software and data, the NLM and the U.S.
  25. *  Government do not and cannot warrant the performance or results that
  26. *  may be obtained by using this software or data. The NLM and the U.S.
  27. *  Government disclaim all warranties, express or implied, including
  28. *  warranties of performance, merchantability or fitness for any particular
  29. *  purpose.
  30. *
  31. *  Please cite the author in any work or product based on this material.
  32. *
  33. * ===========================================================================
  34. *
  35. * Authors:  Aaron Ucko, Wratko Hlavina
  36. *
  37. * File Description:
  38. *   Reader for GFF (including GTF) files.
  39. *
  40. * ===========================================================================
  41. */
  43. #include <ncbi_pch.hpp>
  44. #include <objtools/readers/gff_reader.hpp>
  46. #include <corelib/ncbitime.hpp>
  47. #include <corelib/ncbiutil.hpp>
  49. #include <serial/iterator.hpp>
  51. #include <objects/general/Date.hpp>
  52. #include <objects/seq/Seq_annot.hpp>
  53. #include <objects/seq/Seq_descr.hpp>
  54. #include <objects/seq/Seq_inst.hpp>
  55. #include <objects/seq/Seqdesc.hpp>
  56. #include <objects/seqfeat/Cdregion.hpp>
  57. #include <objects/seqfeat/Gb_qual.hpp>
  58. #include <objects/seqloc/Seq_interval.hpp>
  59. #include <objects/seqloc/Seq_point.hpp>
  60. #include <objects/seqset/Bioseq_set.hpp>
  62. #include <objtools/readers/readfeat.hpp>
  64. #include <algorithm>
  65. #include <ctype.h>
  67. BEGIN_NCBI_SCOPE
  68. BEGIN_SCOPE(objects)
  73. CRef<CSeq_entry> CGFFReader::Read(CNcbiIstream& in, TFlags flags)
  74. {
  75.     x_Reset();
  76.     m_Flags  = flags;
  77.     m_Stream = &in;
  79.     string line;
  80.     while (x_GetNextLine(line)) {
  81.         if (NStr::StartsWith(line, "##")) {
  82.             x_ParseStructuredComment(line);
  83.         } else if (NStr::StartsWith(line, "#")) {
  84.             // regular comment; ignore
  85.         } else {
  86.             CRef<SRecord> record = x_ParseFeatureInterval(line);
  87.             if (record) {
  88.                 string id = x_FeatureID(*record);
  89.                 if (id.empty()) {
  90.                     x_PlaceFeature(*x_ParseRecord(*record), *record);
  91.                 } else {
  92.                     CRef<SRecord>& match = m_DelayedFeats[id];
  93.                     if (match) {
  94.                         x_MergeRecords(*match, *record);
  95.                     } else {
  96.                         match.Reset(record);
  97.                     }
  98.                 }
  99.             }
  100.         }
  101.     }
  103.     ITERATE (TDelayedFeats, it, m_DelayedFeats) {
  104.         x_PlaceFeature(*x_ParseRecord(*it->second), *it->second);
  105.     }
  107.     CRef<CSeq_entry> tse(m_TSE); // need to save before resetting.
  108.     x_Reset();
  109.     return tse;
  110. }
  113. void CGFFReader::x_Warn(const string& message, unsigned int line)
  114. {
  115.     if (line) {
  116.         ERR_POST(Warning << message << " [GFF input, line " << line << ']');
  117.     } else {
  118.         ERR_POST(Warning << message << " [GFF input]");
  119.     }
  120. }
  123. void CGFFReader::x_Reset(void)
  124. {
  125.     m_TSE.Reset(new CSeq_entry);
  126.     m_SeqNameCache.clear();
  127.     m_SeqCache.clear();
  128.     m_DelayedFeats.clear();
  129.     m_DefMol.erase();
  130.     m_LineNumber = 0;
  131. }
  134. bool CGFFReader::x_GetNextLine(string& line)
  135. {
  136.     ++m_LineNumber;
  137.     return NcbiGetlineEOL(*m_Stream, line)  ||  !line.empty();
  138. }
  141. void CGFFReader::x_ParseStructuredComment(const string& line)
  142. {
  143.     vector<string> v;
  144.     NStr::Tokenize(line, "# t", v, NStr::eMergeDelims);
  145.     if (v[0] == "date"  &&  v.size() > 1) {
  146.         x_ParseDateComment(v[1]);
  147.     } else if (v[0] == "type"  &&  v.size() > 1) {
  148.         x_ParseTypeComment(v[1], v.size() > 2 ? v[2] : kEmptyStr);
  149.     }
  150.     // etc.
  151. }
  154. void CGFFReader::x_ParseDateComment(const string& date)
  155. {
  156.     try {
  157.         CRef<CSeqdesc> desc(new CSeqdesc);
  158.         desc->SetUpdate_date().SetToTime(CTime(date, "Y-M-D"),
  159.                                          CDate::ePrecision_day);
  160.         m_TSE->SetSet().SetDescr().Set().push_back(desc);
  161.     } catch (CTimeException& e) {
  162.         x_Warn(string("Bad ISO date: ") + e.what(), x_GetLineNumber());
  163.     }
  164. }
  167. void CGFFReader::x_ParseTypeComment(const string& moltype,
  168.                                     const string& seqname)
  169. {
  170.     if (seqname.empty()) {
  171.         m_DefMol = moltype;
  172.     } else {
  173.         // automatically adds to m_TSE if new
  174.         x_ResolveID(*x_ResolveSeqName(seqname), moltype);
  175.     }
  176. }
  179. CRef<CGFFReader::SRecord>
  180. CGFFReader::x_ParseFeatureInterval(const string& line)
  181. {
  182.     vector<string> v;
  183.     bool           misdelimited = false;
  184.     NStr::Tokenize(line, "t", v, NStr::eNoMergeDelims);
  185.     if (v.size() < 8) {
  186.         NStr::Tokenize(line, " t", v, NStr::eMergeDelims);
  187.         if (v.size() < 8) {
  188.             x_Warn("Skipping line due to insufficient fields",
  189.                    x_GetLineNumber());
  190.             return null;
  191.         } else {
  192.             x_Warn("Bad delimiters (should use tabs)", x_GetLineNumber());
  193.             misdelimited = true;
  194.         }
  195.     } else {
  196.         // XXX - warn about extra fields (if any), but only if they're
  197.         // not comments
  198.         v.resize(9);
  199.     }
  201.     CRef<SRecord> record(x_NewRecord());
  202.     string        accession;
  203.     TSeqPos       from = 0, to = numeric_limits<TSeqPos>::max();
  204.     ENa_strand    strand = eNa_strand_unknown;
  205.     accession      = v[0];
  206.     record->source = v[1];
  207.     record->key    = v[2];
  209.     try {
  210.         from = NStr::StringToUInt(v[3]) - 1;
  211.     } catch (std::exception& e) {
  212.         x_Warn(string("Bad FROM position: ") + e.what(), x_GetLineNumber());
  213.     }
  215.     try {
  216.         to = NStr::StringToUInt(v[4]) - 1;
  217.     } catch (std::exception& e) {
  218.         x_Warn(string("Bad TO position: ") + e.what(), x_GetLineNumber());
  219.     }
  221.     record->score = v[5];
  223.     if (v[6] == "+") {
  224.         strand = eNa_strand_plus;
  225.     } else if (v[6] == "-") {
  226.         strand = eNa_strand_minus;
  227.     } else if (v[6] != ".") {
  228.         x_Warn("Bad strand " + v[6] + " (should be [+-.])", x_GetLineNumber());
  229.     }
  231.     if (v[7] == "0"  ||  v[7] == "1"  ||  v[7] == "2") {
  232.         record->frame = v[7][0] - '0';
  233.     } else if (v[7] == ".") {
  234.         record->frame = -1;
  235.     } else {
  236.         x_Warn("Bad frame " + v[7] + " (should be [012.])", x_GetLineNumber());
  237.         record->frame = -1;
  238.     }
  240.     {{
  241.         SRecord::SSubLoc subloc;
  242.         subloc.accession = accession;
  243.         subloc.strand    = strand;
  244.         subloc.ranges.insert(TSeqRange(from, to));
  246.         record->loc.push_back(subloc);
  247.     }}
  248.    
  249.     string         attr_last_value;
  250.     vector<string> attr_values;
  251.     char           quote_char = 0;
  252.     for (SIZE_TYPE i = 8;  i < v.size();  ++i) {
  253.         string s = v[i] + ' ';
  254.         if ( !misdelimited  &&  i > 8) {
  255.             SIZE_TYPE pos = s.find_first_not_of(" ");
  256.             if (pos != NPOS) {
  257.                 if (s[pos] == '#') {
  258.                     break;
  259.                 } else {
  260.                     x_Warn("Extra non-comment fields", x_GetLineNumber());
  261.                 }
  262.             } else {
  263.                 continue; // just spaces anyway...
  264.             }
  265.         }
  266.         SIZE_TYPE pos = 0;
  267.         while (pos < s.size()) {
  268.             SIZE_TYPE pos2;
  269.             if (quote_char) { // must be inside a value
  270.                 pos2 = s.find_first_of(" '"\", pos);
  271.                 _ASSERT(pos2 != NPOS); // due to trailing space
  272.                 if (s[pos2] == quote_char) {
  273.                     if (attr_values.empty()) {
  274.                         x_Warn("quoted attribute tag " + attr_last_value,
  275.                                x_GetLineNumber());
  276.                     }
  277.                     quote_char = 0;
  278.                     attr_last_value += s.substr(pos, pos2 - pos);
  279.                     try {
  280.                         attr_values.push_back(NStr::ParseEscapes
  281.                                               (attr_last_value));
  282.                     } catch (CStringException& e) {
  283.                         attr_values.push_back(attr_last_value);
  284.                         x_Warn(e.what() + (" in value of " + attr_values[0]),
  285.                                x_GetLineNumber());
  286.                     }
  287.                     attr_last_value.erase();
  288.                 } else if (s[pos2] == '\') {
  289.                     _VERIFY(++pos2 != s.size());
  290.                     attr_last_value += s.substr(pos, pos2 + 1 - pos);
  291.                 } else {
  292.                     attr_last_value += s.substr(pos, pos2 + 1 - pos);
  293.                 }
  294.             } else {
  295.                 pos2 = s.find_first_of(" #;"", pos); // also look for '?
  296.                 _ASSERT(pos2 != NPOS); // due to trailing space
  297.                 if (pos != pos2) {
  298.                     // grab and place the preceding token
  299.                     attr_last_value += s.substr(pos, pos2 - pos);
  300.                     attr_values.push_back(attr_last_value);
  301.                     attr_last_value.erase();
  302.                 }
  304.                 switch (s[pos2]) {
  305.                 case ' ':
  306.                     break;
  308.                 case '#':
  309.                     i = v.size();
  310.                     break;
  312.                 case ';':
  313.                     if (attr_values.empty()) {
  314.                         x_Warn("null attribute", x_GetLineNumber());
  315.                     } else {
  316.                         x_AddAttribute(*record, attr_values);
  317.                         attr_values.clear();
  318.                     }
  319.                     break;
  321.                 // NB: we don't currently search for single quotes.
  322.                 case '"': case ''':
  323.                     quote_char = s[pos2];
  324.                     break;
  326.                 default:
  327.                     _TROUBLE;
  328.                 }
  329.             }
  330.             pos = pos2 + 1;
  331.         }
  332.     }
  334.     if ( !attr_values.empty() ) {
  335.         x_Warn("unterminated attribute " + attr_values[0], x_GetLineNumber());
  336.         x_AddAttribute(*record, attr_values);
  337.     }
  339.     return record;
  340. }
  343. CRef<CSeq_feat> CGFFReader::x_ParseRecord(const SRecord& record)
  344. {
  345.     CRef<CSeq_feat> feat(CFeature_table_reader::CreateSeqFeat
  346.                          (record.key, *x_ResolveLoc(record.loc),
  347.                           CFeature_table_reader::fKeepBadKey));
  348.     if (record.frame >= 0  &&  feat->GetData().IsCdregion()) {
  349.         feat->SetData().SetCdregion().SetFrame
  350.             (static_cast<CCdregion::EFrame>(record.frame + 1));
  351.     }
  352.     ITERATE (SRecord::TAttrs, it, record.attrs) {
  353.         string tag = it->front();
  354.         string value;
  355.         switch (it->size()) {
  356.         case 1:
  357.             break;
  358.         case 2:
  359.             value = (*it)[1];
  360.             break;
  361.         default:
  362.             x_Warn("Ignoring extra fields in value of " + tag, record.line_no);
  363.             value = (*it)[1];
  364.             break;
  365.         }
  366.         if (x_GetFlags() & fGBQuals) {
  367.             if ( !(x_GetFlags() & fNoGTF) ) { // translate
  368.                 if (tag == "transcript_id") {
  369.                     //continue;
  370.                 } else if (tag == "gene_id") {
  371.                     tag = "gene";
  372.                     SIZE_TYPE colon = value.find(':');
  373.                     if (colon != NPOS) {
  374.                         value.erase(0, colon + 1);
  375.                     }
  376.                 } else if (tag == "exon_number") {
  377.                     tag = "number";
  378.                 }
  379.             }
  380.             CFeature_table_reader::AddFeatQual
  381.                 (feat, tag, value, CFeature_table_reader::fKeepBadKey);
  382.         } else { // don't attempt to parse, just treat as imported
  383.             CRef<CGb_qual> qual(new CGb_qual);
  384.             qual->SetQual(tag);
  385.             qual->SetVal(value);
  386.             feat->SetQual().push_back(qual);
  387.         }
  388.     }
  389.     return feat;
  390. }
  393. CRef<CSeq_loc> CGFFReader::x_ResolveLoc(const SRecord::TLoc& loc)
  394. {
  395.     CRef<CSeq_loc> seqloc(new CSeq_loc);
  396.     ITERATE (SRecord::TLoc, it, loc) {
  397.         CRef<CSeq_id> id = x_ResolveSeqName(it->accession);
  398.         ITERATE (set<TSeqRange>, range, it->ranges) {
  399.             CRef<CSeq_loc> segment(new CSeq_loc);
  400.             if (range->GetLength() == 1) {
  401.                 CSeq_point& pnt = segment->SetPnt();
  402.                 pnt.SetId   (*id);
  403.                 pnt.SetPoint(range->GetFrom());
  404.                 if (it->strand != eNa_strand_unknown) {
  405.                     pnt.SetStrand(it->strand);
  406.                 }
  407.             } else {
  408.                 CSeq_interval& si = segment->SetInt();
  409.                 si.SetId  (*id);
  410.                 si.SetFrom(range->GetFrom());
  411.                 si.SetTo  (range->GetTo());
  412.                 if (it->strand != eNa_strand_unknown) {
  413.                     si.SetStrand(it->strand);
  414.                 }
  415.             }
  416.             if (IsReverse(it->strand)) {
  417.                 seqloc->SetMix().Set().push_front(segment);
  418.             } else {
  419.                 seqloc->SetMix().Set().push_back(segment);
  420.             }
  421.         }
  422.     }
  424.     if (seqloc->GetMix().Get().size() == 1) {
  425.         return seqloc->SetMix().Set().front();
  426.     } else {
  427.         return seqloc;
  428.     }
  429. }
  432. void CGFFReader::x_AddAttribute(SRecord& record, vector<string>& attr)
  433. {
  434.     if (x_GetFlags() & fGBQuals) {
  435.         if (attr[0] == "gbkey"  &&  attr.size() == 2) {
  436.             record.key = attr[1];
  437.             return;
  438.         }
  439.     }
  440.     record.attrs.insert(attr);
  441. }
  444. string CGFFReader::x_FeatureID(const SRecord& record)
  445. {
  446.     if (x_GetFlags() & fNoGTF) {
  447.         return kEmptyStr;
  448.     }
  450.     static const vector<string> sc_GeneId(1, "gene_id");
  451.     SRecord::TAttrs::const_iterator gene_it
  452.         = record.attrs.lower_bound(sc_GeneId);
  454.     static const vector<string> sc_TranscriptId(1, "transcript_id");
  455.     SRecord::TAttrs::const_iterator transcript_it
  456.         = record.attrs.lower_bound(sc_TranscriptId);
  458.     static const vector<string> sc_ProteinId(1, "protein_id");
  459.     SRecord::TAttrs::const_iterator protein_it
  460.         = record.attrs.lower_bound(sc_ProteinId);
  462.     // concatenate our IDs from above, if found
  463.     string id;
  464.     if (gene_it != record.attrs.end()  &&
  465.         gene_it->front() == "gene_id") {
  466.         id += (*gene_it)[1];
  467.     }
  469.     if (transcript_it != record.attrs.end()  &&
  470.         transcript_it->front() == "transcript_id") {
  471.         if ( !id.empty() ) {
  472.             id += ' ';
  473.         }
  474.         id += (*transcript_it)[1];
  475.     }
  477.     if (protein_it != record.attrs.end()  &&
  478.         protein_it->front() == "protein_id") {
  479.         if ( !id.empty() ) {
  480.             id += ' ';
  481.         }
  482.         id += (*protein_it)[1];
  483.     }
  485.     // look for db xrefs
  486.     static const vector<string> sc_Dbxref(1, "db_xref");
  487.     SRecord::TAttrs::const_iterator dbxref_it
  488.         = record.attrs.lower_bound(sc_Dbxref);
  489.     for ( ; dbxref_it != record.attrs.end()  &&
  490.             dbxref_it->front() == "db_xref";  ++dbxref_it) {
  491.         if ( !id.empty() ) {
  492.             id += ' ';
  493.         }
  494.         id += (*dbxref_it)[1];
  495.     }
  497.     if ( id.empty() ) {
  498.         return id;
  499.     }
  501.     if (record.key == "start_codon" ||  record.key == "stop_codon") {
  502.         //id += " " + record.key;
  503.         id += "CDS";
  504.     } else if (record.key == "CDS"
  505.                ||  NStr::FindNoCase(record.key, "rna") != NPOS) {
  506.         //id += " " + record.key;
  507.         id += record.key;
  508.     } else { // probably an intron, exon, or single site
  509.         return kEmptyStr;
  510.     }
  511.     _TRACE("id = " << id);
  512.     return id;
  513. }
  516. void CGFFReader::x_MergeRecords(SRecord& dest, const SRecord& src)
  517. {
  518.     // XXX - perform sanity checks and warn on mismatch
  520.     bool merge_overlaps = false;
  521.     if ((src.key == "start_codon"  ||  src.key == "stop_codon")  &&
  522.         dest.key == "CDS") {
  523.         // start_codon and stop_codon features should be contained inside of
  524.         // existing CDS locations
  525.         merge_overlaps = true;
  526.     }
  528.     ITERATE (SRecord::TLoc, slit, src.loc) {
  529.         bool merged = false;
  530.         NON_CONST_ITERATE (SRecord::TLoc, dlit, dest.loc) {
  531.             if (slit->accession != dlit->accession) {
  532.                 if (dest.loc.size() == 1) {
  533.                     x_Warn("Multi-accession feature", src.line_no);
  534.                 }
  535.                 continue;
  536.             } else if (slit->strand != dlit->strand) {
  537.                 if (dest.loc.size() == 1) {
  538.                     x_Warn("Multi-orientation feature", src.line_no);
  539.                 }
  540.                 continue;
  541.             } else {
  542.                 if (merge_overlaps) {
  543.                     ITERATE (set<TSeqRange>, src_iter, slit->ranges) {
  544.                         TSeqRange range(*src_iter);
  545.                         set<TSeqRange>::iterator dst_iter =
  546.                             dlit->ranges.begin();
  547.                         for ( ;  dst_iter != dlit->ranges.end();  ) {
  548.                             if (dst_iter->IntersectingWith(range)) {
  549.                                 range += *dst_iter;
  550.                                 _TRACE("merging overlapping ranges: "
  551.                                        << range.GetFrom() << " - "
  552.                                        << range.GetTo() << " <-> "
  553.                                        << dst_iter->GetFrom() << " - "
  554.                                        << dst_iter->GetTo());
  555.                                 dlit->ranges.erase(dst_iter++);
  556.                             } else {
  557.                                 ++dst_iter;
  558.                             }
  559.                         }
  560.                         dlit->ranges.insert(range);
  561.                     }
  562.                 } else {
  563.                     ITERATE (set<TSeqRange>, set_iter, slit->ranges) {
  564.                         dlit->ranges.insert(*set_iter);
  565.                     }
  566.                 }
  567.                 merged = true;
  568.                 break;
  569.             }
  570.         }
  571.         if ( !merged ) {
  572.             dest.loc.push_back(*slit);
  573.         }
  574.     }
  576.     if (src.key != dest.key) {
  577.         if (dest.key == "CDS"  &&  NStr::EndsWith(src.key, "_codon")
  578.             &&  !(x_GetFlags() & fNoGTF) ) {
  579.             // ok
  580.         } else if (src.key == "CDS" &&  NStr::EndsWith(dest.key, "_codon")
  581.             &&  !(x_GetFlags() & fNoGTF) ) {
  582.             dest.key == "CDS";
  583.         } else {
  584.             x_Warn("Merging features with different keys: " + dest.key
  585.                    + " != " + src.key, src.line_no);
  586.         }
  587.     }
  589.     x_MergeAttributes(dest, src);
  590. }
  593. void CGFFReader::x_MergeAttributes(SRecord& dest, const SRecord& src)
  594. {
  595.     SRecord::TAttrs::iterator dait     = dest.attrs.begin();
  596.     SRecord::TAttrs::iterator dait_end = dest.attrs.end();
  597.     SRecord::TAttrs::iterator dait_tag = dait_end;
  598.     ITERATE (SRecord::TAttrs, sait, src.attrs) {
  599.         const string& tag = sait->front();
  600.         while (dait != dait_end  &&  dait->front() < tag) {
  601.             ++dait;
  602.         }
  603.         if (dait->front() == tag) {
  604.             if (dait_tag == dait_end  ||  dait_tag->front() != tag) {
  605.                 dait_tag = dait;
  606.             }
  607.             while (dait != dait_end  &&  *dait < *sait) {
  608.                 ++dait;
  609.             }
  610.         }
  611.         if (dait != dait_end  &&  *dait == *sait) {
  612.             continue; // identical
  613.         } else if ( !(x_GetFlags() & fNoGTF)  &&  tag == "exon_number") {
  614.             if (dait_tag != dait_end) {
  615.                 while (dait != dait_end  &&  dait->front() == tag) {
  616.                     ++dait;
  617.                 }
  618.                 dest.attrs.erase(dait_tag, dait);
  619.             }
  620.         } else {
  621.             dest.attrs.insert(dait, *sait);
  622.         }
  623.     }
  624. }
  627. void CGFFReader::x_PlaceFeature(CSeq_feat& feat, const SRecord&)
  628. {
  629.     CRef<CBioseq> seq;
  630.     if ( !feat.IsSetProduct() ) {
  631.         for (CTypeConstIterator<CSeq_id> it(feat.GetLocation());  it;  ++it) {
  632.             CRef<CBioseq> seq2 = x_ResolveID(*it, kEmptyStr);
  633.             if ( !seq ) {
  634.                 seq.Reset(seq2);
  635.             } else if ( seq2.NotEmpty()  &&  seq != seq2) {
  636.                 seq.Reset();
  637.                 BREAK(it);
  638.             }
  639.         }
  640.     }
  642.     CBioseq::TAnnot& annots
  643.         = seq ? seq->SetAnnot() : m_TSE->SetSet().SetAnnot();
  644.     NON_CONST_ITERATE (CBioseq::TAnnot, it, annots) {
  645.         if ((*it)->GetData().IsFtable()) {
  646.             (*it)->SetData().SetFtable().push_back(CRef<CSeq_feat>(&feat));
  647.             return;
  648.         }
  649.     }
  650.     CRef<CSeq_annot> annot(new CSeq_annot);
  651.     annot->SetData().SetFtable().push_back(CRef<CSeq_feat>(&feat));
  652.     annots.push_back(annot);
  653. }
  656. CRef<CSeq_id> CGFFReader::x_ResolveSeqName(const string& name)
  657. {
  658.     CRef<CSeq_id>& id = m_SeqNameCache[name];
  659.     if ( !id ) {
  660.         id.Reset(x_ResolveNewSeqName(name));
  661.     }
  662.     if ( !id ) {
  663.         x_Warn("x_ResolveNewSeqName returned null for " + name);
  664.         id.Reset(new CSeq_id(CSeq_id::e_Local, name, name));
  665.     }
  666.     return id;
  667. }
  670. CRef<CSeq_id> CGFFReader::x_ResolveNewSeqName(const string& name)
  671. {
  672.     return CRef<CSeq_id>(new CSeq_id(name));
  673. }
  676. CRef<CBioseq> CGFFReader::x_ResolveID(const CSeq_id& id, const string& mol)
  677. {
  678.     CRef<CBioseq>& seq = m_SeqCache[CConstRef<CSeq_id>(&id)];
  679.     if ( !seq ) {
  680.         seq.Reset(x_ResolveNewID(id, mol));
  681.         // Derived versions of x_ResolveNewID may legimately return null
  682.         // results....
  683.         if (seq) {
  684.             x_PlaceSeq(*seq);
  685.             ITERATE (CBioseq::TId, it, seq->GetId()) {
  686.                 m_SeqCache.insert(make_pair(CConstRef<CSeq_id>(*it), seq));
  687.             }
  688.         }
  689.     }
  690.     return seq;
  691. }
  694. CRef<CBioseq> CGFFReader::x_ResolveNewID(const CSeq_id& id, const string& mol0)
  695. {
  696.     CRef<CBioseq> seq(new CBioseq);
  697.     CRef<CSeq_id> id_copy(new CSeq_id);
  699.     id_copy->Assign(id);
  700.     seq->SetId().push_back(id_copy);
  701.     seq->SetInst().SetRepr(CSeq_inst::eRepr_virtual);
  703.     const string& mol = mol0.empty() ? m_DefMol : mol0;
  704.     if (mol.empty()  ||  mol == "dna") {
  705.         seq->SetInst().SetMol(CSeq_inst::eMol_dna);
  706.     } else if (mol == "rna")  {
  707.         seq->SetInst().SetMol(CSeq_inst::eMol_rna);
  708.     } else if (mol == "protein")  {
  709.         seq->SetInst().SetMol(CSeq_inst::eMol_aa);
  710.     } else {
  711.         x_Warn("unrecognized sequence type " + mol + "; assuming DNA");
  712.         seq->SetInst().SetMol(CSeq_inst::eMol_dna);
  713.     }
  715.     return seq;
  716. }
  719. void CGFFReader::x_PlaceSeq(CBioseq& seq)
  720. {
  721.     bool found = false;
  722.     for (CTypeConstIterator<CBioseq> it(*m_TSE);  it;  ++it) {
  723.         if (&*it == &seq) {
  724.             found = true;
  725.             BREAK(it);
  726.         }
  727.     }
  728.     if ( !found ) {
  729.         CRef<CSeq_entry> se(new CSeq_entry);
  730.         se->SetSeq(seq);
  731.         m_TSE->SetSet().SetSeq_set().push_back(se);
  732.     }
  733. }
  736. END_SCOPE(objects)
  737. END_NCBI_SCOPE
  739. /*
  740. * ===========================================================================
  741. *
  742. * $Log: gff_reader.cpp,v $
  743. * Revision 1000.1  2004/06/01 19:46:20  gouriano
  744. * PRODUCTION: UPGRADED [GCC34_MSVC7] Dev-tree R1.6
  745. *
  746. * Revision 1.6  2004/05/21 21:42:55  gorelenk
  747. * Added PCH ncbi_pch.hpp
  748. *
  749. * Revision 1.5  2004/05/08 12:13:24  dicuccio
  750. * Switched away from using CRangeCollection<> for locations, as this lead to
  751. * inadvertent merging of correctly overlapping exons
  752. *
  753. * Revision 1.4  2004/01/06 17:02:40  dicuccio
  754. * (From Aaron Ucko): Fixed ordering of intervals in a seq-loc-mix on the negative
  755. * strand
  756. *
  757. * Revision 1.3  2003/12/31 12:48:29  dicuccio
  758. * Fixed interpretation of positions - GFF/GTF is 1-based, ASN.1 is 0-based
  759. *
  760. * Revision 1.2  2003/12/04 00:58:24  ucko
  761. * Fix for WorkShop's context-insensitive make_pair.
  762. *
  763. * Revision 1.1  2003/12/03 20:56:36  ucko
  764. * Initial commit.
  765. *
  766. *
  767. * ===========================================================================
  768. */