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

// $Id: BamHeader.cc 3477 2016-03-09 08:26:38Z peter $

/*
  Copyright (C) 2012, 2013, 2014, 2015 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 "yat/utility/split.h"
#include "yat/utility/stl_utility.h"

// we need to include 'sam_header.h' when compiling against libbam
#ifndef YAT_HAVE_HTSLIB
  #if YAT_HAVE_BAM_H
  #include <sam_header.h>
  #elif YAT_HAVE_BAM_BAM_H
  #include <bam/sam_header.h>
  #elif YAT_HAVE_SAMTOOLS_SAM_H
  #include <samtools/sam_header.h>
  #else
  #error cannot end up here
  #endif
#endif

#include <cassert>
#include <cstdlib>
#include <cstring>
#include <string>
#include <sstream>
#include <utility>

namespace theplu {
namespace yat {
namespace omic {

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


  BamHeader::~BamHeader(void)
  {
#if YAT_HAVE_HTSLIB
    bam_hdr_destroy(header_);
#else
    bam_header_destroy(header_);
#endif
  }


  BamHeader::BamHeader(const BamHeader& other)
    : read_group_(other.read_group_), program_group_(other.program_group_)
  {
#ifndef YAT_HAVE_HTSLIB
    using namespace detail; // for bam_hdr_dup in YAT_HAVE_LIBBAM mode
#endif
    header_ = bam_hdr_dup(other.header_);
  }


  const std::string& BamHeader::group(strMap2& table,const std::string& type,
                                      const std::string& id,
                                      const std::string& key) const
  {
    update_group(table, type);
    strMap2::const_iterator line_table = table.find(id);
    if (line_table == table.end()) {
      std::string msg = "@" + type + ": unknown ID: '" + id + "'";
      throw utility::runtime_error(msg);
    }
    strMap::const_iterator it = line_table->second.find(key);
    if (it == line_table->second.end()) {
      std::string msg = "@" + type + ": unknown KEY: '" + key + "'";
      throw utility::runtime_error(msg);
    }
    return it->second;
  }


  void BamHeader::update_group(strMap2& table, const std::string& type) const
  {
    assert(type.size()==2);
    if (table.empty()) { // parse header
      std::istringstream ss(text());
      std::string line;
      while (getline(ss, line)) {
        if (line.substr(1,2)!=type)
          continue;
        // find ID value
        size_t pos = line.find("\tID:", 3);
        size_t end = line.find('\t', pos+4);
        std::string id_val = line.substr(pos+4, end-(pos+4));
        std::vector<std::string> vec;
        utility::split(vec, line, '\t');
        strMap& m = table[id_val];
        for (size_t i=1; i<vec.size(); ++i) {
          m[vec[i].substr(0,2)] = vec[i].substr(3);
        }
      }
    }
  }


  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());
  }


  const std::string& BamHeader::program_group(const std::string& id,
                                              const std::string& key) const
  {
    return group(program_group_, "PG", id, key);
  }


  const std::string& BamHeader::read_group(const std::string& id,
                                           const std::string& key) const
  {
    return group(read_group_, "RG", id, key);
  }


  const std::map<std::string, std::string>&
  BamHeader::read_group(const std::string& id) const
  {
    update_group(read_group_, "RG");
    strMap2::const_iterator line_table = read_group_.find(id);
    if (line_table == read_group_.end()) {
      std::string msg = "@RG: unknown ID: '" + id + "'";
      throw utility::runtime_error(msg);
    }
    return line_table->second;
  }


  BamHeader::read_group_iterator BamHeader::read_group_begin(void) const
  {
    update_group(read_group_, "RG");
    BamHeader::strMap2::const_iterator it(read_group_.begin());
    return utility::pair_first_iterator(it);
  }


  BamHeader::read_group_iterator BamHeader::read_group_end(void) const
  {
    update_group(read_group_, "RG");
    BamHeader::strMap2::const_iterator it(read_group_.end());
    return utility::pair_first_iterator(it);
  }


  void BamHeader::swap(BamHeader& other)
  {
    std::swap(header_, other.header_);
    std::swap(read_group_, other.read_group_);
    std::swap(program_group_, other.program_group_);
  }


  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];
  }


  std::string BamHeader::text(void) const
  {
    assert(header_);
    return std::string(header_->text, header_->l_text);
  }


  void BamHeader::text(const std::string& txt)
  {
    read_group_.clear();
    program_group_.clear();
#ifdef YAT_HAVE_HTSLIB
    // destroy the old
    bam_hdr_destroy(header_);
    // and create header pointer
    header_ = sam_hdr_parse(txt.size(), txt.c_str());
#endif
    header_->l_text = txt.size();
    header_->text = (char*)realloc(header_->text, txt.size()+1);
    memcpy(header_->text, txt.c_str(), txt.size()+1);
#ifndef YAT_HAVE_HTSLIB
    // destroy dict to ensure that things are updated during parsing
    if (header_->dict) {
      sam_header_free(header_->dict);
      header_->dict = 0;
    }
    sam_header_parse(header_);
#endif
  }


  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;
  }


  void swap(BamHeader& lhs, BamHeader& rhs)
  {
    lhs.swap(rhs);
  }


  BamHeader& BamHeader::operator=(const BamHeader& rhs)
  {
    BamHeader tmp(rhs);
    swap(tmp);
    return *this;
  }

#ifndef YAT_HAVE_HTSLIB
  namespace detail {
    bam_header_t * bam_hdr_dup(const bam_header_t* other)
    {
      // taken from htslib 1.2.1 and adjusted to fit old struct bam_header_t
      if (!other)
        return NULL;
      bam_header_t* h = bam_header_init();
      if (!h)
        return NULL;
      // copy stuff
      h->n_targets = other->n_targets;
      h->l_text = other->l_text;
      h->n_text = other->n_text;
      h->dict = NULL;
      h->hash = NULL;
      h->rg2lib = NULL;

      h->text = (char*)calloc(h->l_text + 1, 1);
      memcpy(h->text, other->text, h->l_text);
      h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t));
      h->target_name = (char**)calloc(h->n_targets, sizeof(char*));
      for (int i = 0; i < h->n_targets; ++i) {
        h->target_len[i] = other->target_len[i];
        h->target_name[i] = strdup(other->target_name[i]);
      }
      return h;
    }
  }
#endif

}}}