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

// $Id: BamFile.cc 3408 2015-04-16 00:22:24Z peter $

/*
  Copyright (C) 2012, 2013, 2014 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "BamFile.h"
#include "BamHeader.h"
#include "config_bam.h"

#include "yat/utility/Exception.h"

#include YAT_SAM_HEADER

#include <cassert>
#include <stdexcept>
#include <string>
#include <sstream>

namespace theplu {
namespace yat {
namespace omic {

  InBamFile::InBamFile(void)
    : super_t(), index_(NULL)
  {}


  InBamFile::InBamFile(const std::string& fn)
    : super_t(), index_(NULL)
  {
    open(fn);
  }


  InBamFile::~InBamFile(void)
  {
#if YAT_HAVE_HTSLIB
    hts_idx_destroy(index_);
#else
    bam_index_destroy(index_);
#endif
  }


  const BamHeader& InBamFile::header(void) const
  {
    return header_;
  }


  const InBamFile::index_type* InBamFile::index(void) const
  {
    if (!index_) {
#if YAT_HAVE_HTSLIB
      index_ = hts_idx_load(filename().c_str(), HTS_FMT_BAI);
#else
      index_ = bam_index_load(filename().c_str());
#endif
      if (!index_) {
        std::ostringstream ss;
        ss << "cannot load bam index '" << filename() << ".bai'";
        throw utility::runtime_error(ss.str());
      }
    }
    return index_;
  }


  void InBamFile::open(const std::string& fn)
  {
#if YAT_HAVE_HTSLIB
    open_base(fn, "rb", NULL);
    header_.header_ = sam_hdr_read(sf_);
#else
    open_base(fn, "rb", NULL);
    header_.header_ = detail::bam_hdr_dup(sf_->header);
#endif
    if (header_.header_ == NULL) {
      std::ostringstream ss;
      ss << "failed to read header from '" << filename() << "'";
      throw utility::runtime_error(ss.str());
    }
  }


  bool InBamFile::read(BamRead& read)
  {
#if YAT_HAVE_HTSLIB
    int result_ = sam_read1(sf_, header_.header_, read.bam_);
#else
    int result_ = samread(sf_, read.bam_);
#endif
    if (result_<-1) {
      std::ostringstream ss;
      ss << "InBamFile::read(1): invalid bam file '" << filename() << "'";
      if (result_ == -2)
        ss << " truncated file";
      else
        ss << " error: " << result_;
      throw utility::runtime_error(ss.str());
    }
    return result_>=0;
  }


#if YAT_HAVE_HTSLIB
  bool InBamFile::read(BamRead& read, hts_itr_t* iter)
  {
    int result = sam_itr_next(sf_, iter, read.bam_);
#else
  bool InBamFile::read(BamRead& read, bam_iter_t iter)
  {
    assert(sf_->type & 1); // no random access from sam files
    int result = bam_iter_read(sf_->x.bam, iter, read.bam_);
#endif
    if (result<-1) {
      std::ostringstream ss;
      ss << "read(2): invalid bam file '" << filename() << "'";
      if (result == -2)
        ss << " truncated file";
      else
        ss << " error: " << result;
      throw utility::runtime_error(ss.str());
    }
    return result>=0;
  }


  OutBamFile::OutBamFile(void)
    : super_t() {}


  OutBamFile::OutBamFile(const std::string& fn, const BamHeader& header)
  {
    open(fn, header);
  }


  OutBamFile::OutBamFile(const std::string& fn, const BamHeader& header,
                         unsigned int c)
  {
    open(fn, header, c);
  }


  void OutBamFile::open(const std::string& fn, const BamHeader& h)
  {
#if YAT_HAVE_HTSLIB
    open_base(fn, "wb", NULL);
    sam_hdr_write(sf_, h.header_);
#else
    open_base(fn, "wb", h.header_);
#endif
  }


  void OutBamFile::open(const std::string& fn, const BamHeader& h,
                        unsigned int compression)
  {
    std::string mode("wb");
    if (compression > 9) {
      std::stringstream oss;
      oss << "OutBamFile::open( " << fn << ", <header>, " << compression
          << "): incorrect compression level\n";
      throw std::invalid_argument(oss.str());
    }
#if YAT_HAVE_HTSLIB
    if (compression)
      mode.push_back('0'+compression);
    else
      mode = "wu";
    open_base(fn, mode, NULL);
    sam_hdr_write(sf_, h.header_);
#else
    mode.push_back('0'+compression);
    open_base(fn, mode, h.header_);
#endif
  }


  void OutBamFile::write(const BamRead& read)
  {
#if YAT_HAVE_HTSLIB
    assert(read.bam_);
    if (is_open() && bam_write1(sf_->fp.bgzf, read.bam_)>0)
      return;
#else
    if (samwrite(sf_, read.bam_)>0)
      return;
#endif
    throw OutBamFile::error(read);
  }


  // OutBamFile::error

  OutBamFile::error::error(const BamRead& r)
    : IO_error(""), read_(r) {}


  OutBamFile::error::~error(void) throw() {}


  const BamRead& OutBamFile::error::read(void) const
  { return read_; }

}}}