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

Last change on this file since 3408 was 3408, checked in by Peter, 4 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.7 KB
RevLine 
[2883]1// $Id: BamFile.cc 3408 2015-04-16 00:22:24Z peter $
2
[2993]3/*
[3160]4  Copyright (C) 2012, 2013, 2014 Peter Johansson
[2993]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
[2883]22#include <config.h>
23
24#include "BamFile.h"
25#include "BamHeader.h"
[2943]26#include "config_bam.h"
[2883]27
[2897]28#include "yat/utility/Exception.h"
29
[3350]30#include YAT_SAM_HEADER
[2883]31
[3048]32#include <cassert>
[2883]33#include <stdexcept>
34#include <string>
35#include <sstream>
36
37namespace theplu {
38namespace yat {
39namespace omic {
40
41  InBamFile::InBamFile(void)
[2895]42    : super_t(), index_(NULL)
[2883]43  {}
44
45
46  InBamFile::InBamFile(const std::string& fn)
[2895]47    : super_t(), index_(NULL)
[2883]48  {
49    open(fn);
50  }
51
52
53  InBamFile::~InBamFile(void)
54  {
[3353]55#if YAT_HAVE_HTSLIB
56    hts_idx_destroy(index_);
57#else
[2895]58    bam_index_destroy(index_);
[3353]59#endif
[2883]60  }
61
62
63  const BamHeader& InBamFile::header(void) const
64  {
65    return header_;
66  }
67
68
[3353]69  const InBamFile::index_type* InBamFile::index(void) const
[2883]70  {
71    if (!index_) {
[3353]72#if YAT_HAVE_HTSLIB
73      index_ = hts_idx_load(filename().c_str(), HTS_FMT_BAI);
74#else
[2883]75      index_ = bam_index_load(filename().c_str());
[3353]76#endif
[2883]77      if (!index_) {
78        std::ostringstream ss;
79        ss << "cannot load bam index '" << filename() << ".bai'";
[2897]80        throw utility::runtime_error(ss.str());
[2883]81      }
82    }
83    return index_;
84  }
85
86
87  void InBamFile::open(const std::string& fn)
88  {
[3353]89#if YAT_HAVE_HTSLIB
[3357]90    open_base(fn, "rb", NULL);
91    header_.header_ = sam_hdr_read(sf_);
[3353]92#else
[2883]93    open_base(fn, "rb", NULL);
[3408]94    header_.header_ = detail::bam_hdr_dup(sf_->header);
[3353]95#endif
[3079]96    if (header_.header_ == NULL) {
97      std::ostringstream ss;
98      ss << "failed to read header from '" << filename() << "'";
99      throw utility::runtime_error(ss.str());
100    }
[2883]101  }
102
103
104  bool InBamFile::read(BamRead& read)
105  {
[3353]106#if YAT_HAVE_HTSLIB
107    int result_ = sam_read1(sf_, header_.header_, read.bam_);
108#else
[2883]109    int result_ = samread(sf_, read.bam_);
[3353]110#endif
[2883]111    if (result_<-1) {
112      std::ostringstream ss;
[3080]113      ss << "InBamFile::read(1): invalid bam file '" << filename() << "'";
[2883]114      if (result_ == -2)
115        ss << " truncated file";
116      else
117        ss << " error: " << result_;
[2897]118      throw utility::runtime_error(ss.str());
[2883]119    }
120    return result_>=0;
121  }
122
123
[3353]124#if YAT_HAVE_HTSLIB
125  bool InBamFile::read(BamRead& read, hts_itr_t* iter)
126  {
[3358]127    int result = sam_itr_next(sf_, iter, read.bam_);
[3353]128#else
[2883]129  bool InBamFile::read(BamRead& read, bam_iter_t iter)
130  {
131    assert(sf_->type & 1); // no random access from sam files
[3358]132    int result = bam_iter_read(sf_->x.bam, iter, read.bam_);
[3353]133#endif
[3358]134    if (result<-1) {
[2883]135      std::ostringstream ss;
136      ss << "read(2): invalid bam file '" << filename() << "'";
[3358]137      if (result == -2)
[2883]138        ss << " truncated file";
139      else
[3358]140        ss << " error: " << result;
[2897]141      throw utility::runtime_error(ss.str());
[2883]142    }
[3358]143    return result>=0;
[2883]144  }
145
146
147  OutBamFile::OutBamFile(void)
148    : super_t() {}
149
150
151  OutBamFile::OutBamFile(const std::string& fn, const BamHeader& header)
152  {
153    open(fn, header);
154  }
155
156
[3160]157  OutBamFile::OutBamFile(const std::string& fn, const BamHeader& header,
158                         unsigned int c)
159  {
160    open(fn, header, c);
161  }
162
163
[2883]164  void OutBamFile::open(const std::string& fn, const BamHeader& h)
165  {
[3353]166#if YAT_HAVE_HTSLIB
[3357]167    open_base(fn, "wb", NULL);
168    sam_hdr_write(sf_, h.header_);
[3353]169#else
[2883]170    open_base(fn, "wb", h.header_);
[3353]171#endif
[2883]172  }
173
174
[3160]175  void OutBamFile::open(const std::string& fn, const BamHeader& h,
176                        unsigned int compression)
177  {
178    std::string mode("wb");
179    if (compression > 9) {
180      std::stringstream oss;
181      oss << "OutBamFile::open( " << fn << ", <header>, " << compression
182          << "): incorrect compression level\n";
183      throw std::invalid_argument(oss.str());
184    }
[3353]185#if YAT_HAVE_HTSLIB
[3357]186    if (compression)
187      mode.push_back('0'+compression);
188    else
189      mode = "wu";
190    open_base(fn, mode, NULL);
191    sam_hdr_write(sf_, h.header_);
[3353]192#else
[3160]193    mode.push_back('0'+compression);
194    open_base(fn, mode, h.header_);
[3353]195#endif
[3160]196  }
197
198
[2883]199  void OutBamFile::write(const BamRead& read)
200  {
[3353]201#if YAT_HAVE_HTSLIB
[3359]202    assert(read.bam_);
203    if (is_open() && bam_write1(sf_->fp.bgzf, read.bam_)>0)
204      return;
[3353]205#else
[3359]206    if (samwrite(sf_, read.bam_)>0)
207      return;
[3353]208#endif
[3359]209    throw OutBamFile::error(read);
[2883]210  }
211
[3166]212
213  // OutBamFile::error
214
215  OutBamFile::error::error(const BamRead& r)
216    : IO_error(""), read_(r) {}
217
218
219  OutBamFile::error::~error(void) throw() {}
220
221
222  const BamRead& OutBamFile::error::read(void) const
223  { return read_; }
224
[2883]225}}}
Note: See TracBrowser for help on using the repository browser.