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

Last change on this file since 3359 was 3359, checked in by Peter, 5 years ago

fix issues in OutBamFile::write. refs #794

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 4.7 KB
Line 
1// $Id: BamFile.cc 3359 2014-11-23 12:05:13Z 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 "BamFile.h"
25#include "BamHeader.h"
26#include "config_bam.h"
27
28#include "yat/utility/Exception.h"
29
30#include YAT_SAM_HEADER
31
32#include <cassert>
33#include <stdexcept>
34#include <string>
35#include <sstream>
36
37namespace theplu {
38namespace yat {
39namespace omic {
40
41  InBamFile::InBamFile(void)
42    : super_t(), index_(NULL)
43  {}
44
45
46  InBamFile::InBamFile(const std::string& fn)
47    : super_t(), index_(NULL)
48  {
49    open(fn);
50  }
51
52
53  InBamFile::~InBamFile(void)
54  {
55#if YAT_HAVE_HTSLIB
56    hts_idx_destroy(index_);
57#else
58    bam_index_destroy(index_);
59#endif
60  }
61
62
63  const BamHeader& InBamFile::header(void) const
64  {
65    return header_;
66  }
67
68
69  const InBamFile::index_type* InBamFile::index(void) const
70  {
71    if (!index_) {
72#if YAT_HAVE_HTSLIB
73      index_ = hts_idx_load(filename().c_str(), HTS_FMT_BAI);
74#else
75      index_ = bam_index_load(filename().c_str());
76#endif
77      if (!index_) {
78        std::ostringstream ss;
79        ss << "cannot load bam index '" << filename() << ".bai'";
80        throw utility::runtime_error(ss.str());
81      }
82    }
83    return index_;
84  }
85
86
87  void InBamFile::open(const std::string& fn)
88  {
89#if YAT_HAVE_HTSLIB
90    open_base(fn, "rb", NULL);
91    header_.header_ = sam_hdr_read(sf_);
92#else
93    open_base(fn, "rb", NULL);
94    header_.header_ = sf_->header;
95#endif
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    }
101  }
102
103
104  bool InBamFile::read(BamRead& read)
105  {
106#if YAT_HAVE_HTSLIB
107    int result_ = sam_read1(sf_, header_.header_, read.bam_);
108#else
109    int result_ = samread(sf_, read.bam_);
110#endif
111    if (result_<-1) {
112      std::ostringstream ss;
113      ss << "InBamFile::read(1): invalid bam file '" << filename() << "'";
114      if (result_ == -2)
115        ss << " truncated file";
116      else
117        ss << " error: " << result_;
118      throw utility::runtime_error(ss.str());
119    }
120    return result_>=0;
121  }
122
123
124#if YAT_HAVE_HTSLIB
125  bool InBamFile::read(BamRead& read, hts_itr_t* iter)
126  {
127    int result = sam_itr_next(sf_, iter, read.bam_);
128#else
129  bool InBamFile::read(BamRead& read, bam_iter_t iter)
130  {
131    assert(sf_->type & 1); // no random access from sam files
132    int result = bam_iter_read(sf_->x.bam, iter, read.bam_);
133#endif
134    if (result<-1) {
135      std::ostringstream ss;
136      ss << "read(2): invalid bam file '" << filename() << "'";
137      if (result == -2)
138        ss << " truncated file";
139      else
140        ss << " error: " << result;
141      throw utility::runtime_error(ss.str());
142    }
143    return result>=0;
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
157  OutBamFile::OutBamFile(const std::string& fn, const BamHeader& header,
158                         unsigned int c)
159  {
160    open(fn, header, c);
161  }
162
163
164  void OutBamFile::open(const std::string& fn, const BamHeader& h)
165  {
166#if YAT_HAVE_HTSLIB
167    open_base(fn, "wb", NULL);
168    sam_hdr_write(sf_, h.header_);
169#else
170    open_base(fn, "wb", h.header_);
171#endif
172  }
173
174
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    }
185#if YAT_HAVE_HTSLIB
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_);
192#else
193    mode.push_back('0'+compression);
194    open_base(fn, mode, h.header_);
195#endif
196  }
197
198
199  void OutBamFile::write(const BamRead& read)
200  {
201#if YAT_HAVE_HTSLIB
202    assert(read.bam_);
203    if (is_open() && bam_write1(sf_->fp.bgzf, read.bam_)>0)
204      return;
205#else
206    if (samwrite(sf_, read.bam_)>0)
207      return;
208#endif
209    throw OutBamFile::error(read);
210  }
211
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
225}}}
Note: See TracBrowser for help on using the repository browser.