source: trunk/yat/omic/BamHeader.cc @ 3409

Last change on this file since 3409 was 3409, checked in by Peter, 8 years ago

refs #817. Fix bug that struct was not updated in text(1) function (in libbam mode); add test to catch bug.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 5.4 KB
Line 
1// $Id: BamHeader.cc 3409 2015-04-17 01:07:04Z peter $
2
3/*
4  Copyright (C) 2012, 2013, 2014 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <config.h>
23
24#include "BamHeader.h"
25
26#include "yat/utility/Exception.h"
27
28// we need to include 'sam_header.h' when compiling against libbam
29#ifndef YAT_HAVE_HTSLIB
30  #if YAT_HAVE_SAM_H
31  #include <sam_header.h>
32  #elif YAT_HAVE_BAM_BAM_H
33  #include <bam/sam_header.h>
34  #elif YAT_HAVE_SAMTOOLS_SAM_H
35  #include <samtools/sam_header.h>
36  #else
37  #error cannot end up here
38  #endif
39#endif
40
41#include <cassert>
42#include <cstdlib>
43#include <cstring>
44#include <string>
45#include <sstream>
46#include <utility>
47
48namespace theplu {
49namespace yat {
50namespace omic {
51
52  BamHeader::BamHeader(void)
53    : header_(NULL)
54  {
55  }
56
57
58  BamHeader::~BamHeader(void)
59  {
60#if YAT_HAVE_HTSLIB
61    bam_hdr_destroy(header_);
62#else
63    bam_header_destroy(header_);
64#endif
65  }
66
67
68  BamHeader::BamHeader(const BamHeader& other)
69  {
70#ifndef YAT_HAVE_HTSLIB
71    using namespace detail; // for bam_hdr_dup in YAT_HAVE_LIBBAM mode
72#endif
73    header_ = bam_hdr_dup(other.header_);
74  }
75
76
77  void BamHeader::parse_region(const std::string& reg, int& id, int& begin,
78                               int& end) const
79  {
80    assert(header_);
81#if YAT_HAVE_HTSLIB
82    const char* b = reg.c_str();
83    const char* e = hts_parse_reg(b, &begin, &end);
84    // If begin > end suggests something is wrong and that is also how
85    // bam_parse_region (below in samtools impl) used to detect error
86    // and return non-zero.
87    if (begin<=end) {
88      std::string chr = reg.substr(0, e-b);
89      try {
90        id = tid(chr);
91        return;
92      }
93      catch (utility::runtime_error& e) {
94        ; // throw below instead
95      }
96    }
97#else
98    if (!bam_parse_region(header_, reg.c_str(), &id, &begin, &end))
99      return;
100#endif
101    std::ostringstream ss;
102    ss << "invalid region: '" << reg << "'";
103    throw utility::runtime_error(ss.str());
104  }
105
106
107  void BamHeader::swap(BamHeader& other)
108  {
109    std::swap(header_, other.header_);
110  }
111
112
113  uint32_t BamHeader::target_length(size_t tid) const
114  {
115    assert(tid < static_cast<size_t>(n_targets()));
116    return header_->target_len[tid];
117  }
118
119
120  const char* BamHeader::target_name(size_t tid) const
121  {
122    assert(tid < static_cast<size_t>(n_targets()));
123    return header_->target_name[tid];
124  }
125
126
127  std::string BamHeader::text(void) const
128  {
129    assert(header_);
130    return std::string(header_->text, header_->l_text);
131  }
132
133
134  void BamHeader::text(const std::string& txt)
135  {
136#ifdef YAT_HAVE_HTSLIB
137    // destroy the old
138    bam_hdr_destroy(header_);
139    // and create header pointer
140    header_ = sam_hdr_parse(txt.size(), txt.c_str());
141#endif
142    header_->l_text = txt.size();
143    header_->text = (char*)realloc(header_->text, txt.size()+1);
144    memcpy(header_->text, txt.c_str(), txt.size()+1);
145#ifndef YAT_HAVE_HTSLIB
146    // destroy dict to ensure that things are updated during parsing
147    if (header_->dict) {
148      sam_header_free(header_->dict);
149      header_->dict = 0;
150    }
151    sam_header_parse(header_);
152#endif
153  }
154
155
156  int32_t BamHeader::tid(const std::string& name) const
157  {
158#if YAT_HAVE_HTSLIB
159    int res = bam_name2id(header_, name.c_str());
160    if (res==-1) {
161      std::string msg("invalid segment name:");
162      msg += name;
163      throw utility::runtime_error(msg);
164    }
165    return res;
166#else
167    if (!header_->hash) {
168      // We would like to call something like bam_init_header_hash,
169      // but unfortunately that function is not defined in any header
170      // (only in two source files). Instead we call bam_parse_region
171      // because a documented side-effect of that function is that
172      // hash get populated.
173      int tid, b, e;
174      bam_parse_region(header_, name.c_str(), &tid, &b, &e);
175      assert(header_->hash);
176      assert(tid == bam_get_tid(header_, name.c_str()));
177      return tid;
178    }
179    return bam_get_tid(header_, name.c_str());
180#endif
181  }
182
183
184  int32_t BamHeader::n_targets(void) const
185  {
186    assert(header_);
187    return header_->n_targets;
188  }
189
190
191  void swap(BamHeader& lhs, BamHeader& rhs)
192  {
193    lhs.swap(rhs);
194  }
195
196
197  BamHeader& BamHeader::operator=(const BamHeader& rhs)
198  {
199    BamHeader tmp(rhs);
200    swap(tmp);
201    return *this;
202  }
203
204#ifndef YAT_HAVE_HTSLIB
205  namespace detail {
206    bam_header_t * bam_hdr_dup(const bam_header_t* other)
207    {
208      // taken from htslib 1.2.1 and adjusted to fit old struct bam_header_t
209      if (!other)
210        return NULL;
211      bam_header_t* h = bam_header_init();
212      if (!h)
213        return NULL;
214      // copy stuff
215      h->n_targets = other->n_targets;
216      h->l_text = other->l_text;
217      h->n_text = other->n_text;
218      h->dict = NULL;
219      h->hash = NULL;
220      h->rg2lib = NULL;
221
222      h->text = (char*)calloc(h->l_text + 1, 1);
223      memcpy(h->text, other->text, h->l_text);
224      h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t));
225      h->target_name = (char**)calloc(h->n_targets, sizeof(char*));
226      for (int i = 0; i < h->n_targets; ++i) {
227        h->target_len[i] = other->target_len[i];
228        h->target_name[i] = strdup(other->target_name[i]);
229      }
230      return h;
231    }
232  }
233#endif
234
235}}}
Note: See TracBrowser for help on using the repository browser.