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

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

Implement functions to get/set free text in BamHeader?. As setting the
header text modifies underlying data, copy and assignment are now
implemented as hard copy (rather than the implicit pointer copy as
before). It also means that the bam_hdr_t* is no longer owned by
InBamFile?, but is now owned by BamHeader? and destroyed in its
destructor.

refs #817

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.9 KB
Line 
1// $Id: BamHeader.cc 3408 2015-04-16 00:22:24Z 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#include <cassert>
29#include <cstdlib>
30#include <cstring>
31#include <string>
32#include <sstream>
33#include <utility>
34
35namespace theplu {
36namespace yat {
37namespace omic {
38
39  BamHeader::BamHeader(void)
40    : header_(NULL)
41  {
42  }
43
44
45  BamHeader::~BamHeader(void)
46  {
47#if YAT_HAVE_HTSLIB
48    bam_hdr_destroy(header_);
49#else
50    bam_header_destroy(header_);
51#endif
52  }
53
54
55  BamHeader::BamHeader(const BamHeader& other)
56  {
57#ifndef YAT_HAVE_HTSLIB
58    using namespace detail; // for bam_hdr_dup in YAT_HAVE_LIBBAM mode
59#endif
60    header_ = bam_hdr_dup(other.header_);
61  }
62
63
64  void BamHeader::parse_region(const std::string& reg, int& id, int& begin,
65                               int& end) const
66  {
67    assert(header_);
68#if YAT_HAVE_HTSLIB
69    const char* b = reg.c_str();
70    const char* e = hts_parse_reg(b, &begin, &end);
71    // If begin > end suggests something is wrong and that is also how
72    // bam_parse_region (below in samtools impl) used to detect error
73    // and return non-zero.
74    if (begin<=end) {
75      std::string chr = reg.substr(0, e-b);
76      try {
77        id = tid(chr);
78        return;
79      }
80      catch (utility::runtime_error& e) {
81        ; // throw below instead
82      }
83    }
84#else
85    if (!bam_parse_region(header_, reg.c_str(), &id, &begin, &end))
86      return;
87#endif
88    std::ostringstream ss;
89    ss << "invalid region: '" << reg << "'";
90    throw utility::runtime_error(ss.str());
91  }
92
93
94  void BamHeader::swap(BamHeader& other)
95  {
96    std::swap(header_, other.header_);
97  }
98
99
100  uint32_t BamHeader::target_length(size_t tid) const
101  {
102    assert(tid < static_cast<size_t>(n_targets()));
103    return header_->target_len[tid];
104  }
105
106
107  const char* BamHeader::target_name(size_t tid) const
108  {
109    assert(tid < static_cast<size_t>(n_targets()));
110    return header_->target_name[tid];
111  }
112
113
114  std::string BamHeader::text(void) const
115  {
116    assert(header_);
117    return std::string(header_->text, header_->l_text);
118  }
119
120  void BamHeader::text(const std::string& txt)
121  {
122#ifdef YAT_HAVE_HTSLIB
123    // destroy the old
124    bam_hdr_destroy(header_);
125    // and create header pointer
126    header_ = sam_hdr_parse(txt.size(), txt.c_str());
127#endif
128    header_->l_text = txt.size();
129    header_->text = (char*)realloc(header_->text, txt.size()+1);
130    memcpy(header_->text, txt.c_str(), txt.size()+1);
131#ifndef YAT_HAVE_HTSLIB
132    sam_header_parse(header_);
133#endif
134  }
135
136
137  int32_t BamHeader::tid(const std::string& name) const
138  {
139#if YAT_HAVE_HTSLIB
140    int res = bam_name2id(header_, name.c_str());
141    if (res==-1) {
142      std::string msg("invalid segment name:");
143      msg += name;
144      throw utility::runtime_error(msg);
145    }
146    return res;
147#else
148    if (!header_->hash) {
149      // We would like to call something like bam_init_header_hash,
150      // but unfortunately that function is not defined in any header
151      // (only in two source files). Instead we call bam_parse_region
152      // because a documented side-effect of that function is that
153      // hash get populated.
154      int tid, b, e;
155      bam_parse_region(header_, name.c_str(), &tid, &b, &e);
156      assert(header_->hash);
157      assert(tid == bam_get_tid(header_, name.c_str()));
158      return tid;
159    }
160    return bam_get_tid(header_, name.c_str());
161#endif
162  }
163
164
165  int32_t BamHeader::n_targets(void) const
166  {
167    assert(header_);
168    return header_->n_targets;
169  }
170
171
172  void swap(BamHeader& lhs, BamHeader& rhs)
173  {
174    lhs.swap(rhs);
175  }
176
177
178  BamHeader& BamHeader::operator=(const BamHeader& rhs)
179  {
180    BamHeader tmp(rhs);
181    swap(tmp);
182    return *this;
183  }
184
185#ifndef YAT_HAVE_HTSLIB
186  namespace detail {
187    bam_header_t * bam_hdr_dup(const bam_header_t* other)
188    {
189      // taken from htslib 1.2.1 and adjusted to fit old struct bam_header_t
190      if (!other)
191        return NULL;
192      bam_header_t* h = bam_header_init();
193      if (!h)
194        return NULL;
195      // copy stuff
196      h->n_targets = other->n_targets;
197      h->l_text = other->l_text;
198      h->n_text = other->n_text;
199      h->dict = NULL;
200      h->hash = NULL;
201      h->rg2lib = NULL;
202
203      h->text = (char*)calloc(h->l_text + 1, 1);
204      memcpy(h->text, other->text, h->l_text);
205      h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t));
206      h->target_name = (char**)calloc(h->n_targets, sizeof(char*));
207      for (int i = 0; i < h->n_targets; ++i) {
208        h->target_len[i] = other->target_len[i];
209        h->target_name[i] = strdup(other->target_name[i]);
210      }
211      return h;
212    }
213  }
214#endif
215
216}}}
Note: See TracBrowser for help on using the repository browser.