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

// $Id: BamHeader.cc 3360 2014-11-24 01:22:23Z 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 "BamHeader.h"

#include "yat/utility/Exception.h"

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

namespace theplu {
namespace yat {
namespace omic {

  BamHeader::BamHeader(void)
    : header_(NULL)
  {
  }


  void BamHeader::parse_region(const std::string& reg, int& id, int& begin,
                               int& end) const
  {
    assert(header_);
#if YAT_HAVE_HTSLIB
    const char* b = reg.c_str();
    const char* e = hts_parse_reg(b, &begin, &end);
    // If begin > end suggests something is wrong and that is also how
    // bam_parse_region (below in samtools impl) used to detect error
    // and return non-zero.
    if (begin<=end) {
      std::string chr = reg.substr(0, e-b);
      try {
        id = tid(chr);
        return;
      }
      catch (utility::runtime_error& e) {
        ; // throw below instead
      }
    }
#else
    if (!bam_parse_region(header_, reg.c_str(), &id, &begin, &end))
      return;
#endif
    std::ostringstream ss;
    ss << "invalid region: '" << reg << "'";
    throw utility::runtime_error(ss.str());
  }


  uint32_t BamHeader::target_length(size_t tid) const
  {
    assert(tid < static_cast<size_t>(n_targets()));
    return header_->target_len[tid];
  }


  const char* BamHeader::target_name(size_t tid) const
  {
    assert(tid < static_cast<size_t>(n_targets()));
    return header_->target_name[tid];
  }


  int32_t BamHeader::tid(const std::string& name) const
  {
#if YAT_HAVE_HTSLIB
    int res = bam_name2id(header_, name.c_str());
    if (res==-1) {
      std::string msg("invalid segment name:");
      msg += name;
      throw utility::runtime_error(msg);
    }
    return res;
#else
    if (!header_->hash) {
      // We would like to call something like bam_init_header_hash,
      // but unfortunately that function is not defined in any header
      // (only in two source files). Instead we call bam_parse_region
      // because a documented side-effect of that function is that
      // hash get populated.
      int tid, b, e;
      bam_parse_region(header_, name.c_str(), &tid, &b, &e);
      assert(header_->hash);
      assert(tid == bam_get_tid(header_, name.c_str()));
      return tid;
    }
    return bam_get_tid(header_, name.c_str());
#endif
  }


  int32_t BamHeader::n_targets(void) const
  {
    assert(header_);
    return header_->n_targets;
  }

}}}