source: trunk/yat/omic/BamRead.cc @ 3031

// $Id: BamRead.cc 3031 2013-04-27 11:24:13Z peter $

/*
  Copyright (C) 2012, 2013 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 "BamRead.h"
#include "config_bam.h"

#include YAT_BAM_HEADER

#include "yat/utility/Exception.h"

#include <algorithm>
#include <cassert>
#include <sstream>

// bam1_seq_seti is not defined in 0.1.18 (or older), so
// backport. Code is taken bam.h in version 0.1.19-dev.
#ifndef bam1_seq_seti
#define bam1_seq_seti(s, i, c) ( (s)[(i)>>1] = \
  ((s)[(i)>>1] & 0xf<<(((i)&1)<<2)) | (c)<<((~(i)&1)<<2) )
#endif

namespace theplu {
namespace yat {
namespace omic {

  BamRead::BamRead(void)
    : bam_(bam_init1())
  {}


  BamRead::BamRead(const BamRead& other)
    : bam_(bam_dup1(other.bam_))
  {
  }


  BamRead::~BamRead(void)
  {
    bam_destroy1(bam_);
  }


  const BamRead& BamRead::operator=(const BamRead& rhs)
  {
    if (bam_!=rhs.bam_)
      bam_copy1(bam_, rhs.bam_);
    return *this;
  }


  const uint8_t* BamRead::aux(void) const
  {
    assert(bam_);
    return bam1_aux(bam_);
  }


  const uint8_t* BamRead::aux(const char tag[2]) const
  {
    assert(bam_);
    return bam_aux_get(bam_, tag);
  }


  void BamRead::aux_append(const char tag[2], char type,int len,uint8_t* data)
  {
    assert(bam_);
    bam_aux_append(bam_, tag, type, len, data);
  }


  void BamRead::aux_del(const char tag[2])
  {
    assert(bam_);
    uint8_t* s = bam_aux_get(bam_, tag);
    if (!s) {
      std::stringstream msg("BamRead::aux_del: absent tag: ");
      msg << "'" << tag << "'";
      throw utility::runtime_error(msg.str());
    }
    bam_aux_del(bam_, s);
  }


  int BamRead::aux_size(void) const
  {
    assert(bam_);
    return bam_->l_aux;
  }


  const uint32_t* BamRead::cigar(void) const
  {
    return bam1_cigar(bam_);
  }


  uint32_t BamRead::cigar(size_t i) const
  {
    assert(i<core().n_cigar);
    return bam1_cigar(bam_)[i];
  }


  uint32_t BamRead::cigar_op(size_t i) const
  {
    assert(i<core().n_cigar);
    return bam_cigar_op(cigar(i));
  }


  uint32_t BamRead::cigar_oplen(size_t i) const
  {
    assert(i<core().n_cigar);
    return bam_cigar_oplen(cigar(i));
  }


  std::string BamRead::cigar_str(void) const
  {
    std::ostringstream os;
    std::string str = BAM_CIGAR_STR;
    for (uint32_t i=0; i<core().n_cigar; ++i) {
      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
         << str[bam_cigar_op(bam1_cigar(bam_)[i])];
    }
    return os.str();
  }


  void BamRead::cigar(const std::vector<uint32_t>& c)
  {
    // special case when cigar size is not changed
    if (c.size() == core().n_cigar) {
      if (c.size())
        memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
      return;
    }

    int m_data = bam_->m_data;
    int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
    // double capacity if needed
    if (m_data < data_len) {
      m_data = data_len;
      kroundup32(m_data);
    }
    uint8_t* data = (uint8_t*)calloc(m_data, 1);
    /*
      data comes as qname, cigar, rest
     */
    // first copy qname
    memcpy(data, bam_->data, core().l_qname);
    // copy new cigar
    if (c.size())
      // each cigar element is 4 bytes
      memcpy(data + core().l_qname, &c[0], 4*c.size());
    // copy remaining data
    memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
           data_len - 4*c.size() - core().l_qname);
    bam_->m_data = m_data;
    bam_->data_len = data_len;
    core().n_cigar = c.size();
    free(bam_->data);
    bam_->data = data;
  }


  const bam1_core_t& BamRead::core(void) const
  { return bam_->core; }


  bam1_core_t& BamRead::core(void)
  { return bam_->core; }


  const char* BamRead::name(void) const
  { return bam1_qname(bam_); }


  uint8_t BamRead::qual(size_t i) const
  { return *(bam1_qual(bam_)+i); }


  void BamRead::qual(size_t i, uint8_t q)
  { *(bam1_qual(bam_)+i)=q; }


  void BamRead::reserve(int n)
  {
    assert(bam_);
    if (bam_->m_data >= n)
      return;
    // at least double capacity
    bam_->m_data = std::max(bam_->m_data<<1, n);
    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
    if (!bam_->data) {
      throw utility::runtime_error("BamRead::reserve failed");
    }
  }


  std::string BamRead::sequence(void) const
  {
    std::string result(sequence_length(), ' ');
#ifdef HAVE_BAM_NT16_REV_TABLE
    for (size_t i=0; i<result.size(); ++i)
      result[i] = bam_nt16_rev_table[sequence(i)];
#else
    std::string table = "=ACMGRSVTWYHKDBN";
    for (size_t i=0; i<result.size(); ++i)
      result[i] = table[sequence(i)];
#endif
    return result;
  }


  void BamRead::sequence(const std::string& seq,
                         const std::vector<uint8_t>& qual)
  {
    assert(seq.size()==qual.size());
    // data consist of
    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->

    uint8_t* new_aux = bam1_seq(bam_) + (seq.size()+1)/2 + qual.size();

    int len = new_aux + bam_->l_aux - bam_->data;
    reserve(len);
    // move aux to its new location
    memmove(static_cast<void*>(new_aux),
            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
    // copy sequence
    for (size_t i=0; i<seq.size(); ++i)
      sequence(i, bam_nt16_table[seq[i]]);
    core().l_qseq = seq.size();
    // copy quality
    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
      memcpy(static_cast<void*>(bam1_qual(bam_)),
             static_cast<const void*>(&qual[0]), qual.size());
  }


  void BamRead::sequence(size_t i, uint8_t seq)
  {
    assert(bam_);
    assert(seq < 16);
    bam1_seq_seti(bam1_seq(bam_), i, seq);
  }


  uint8_t BamRead::sequence(size_t index) const
  {
    assert(bam_);
    return bam1_seqi(bam1_seq(bam_),index);
  }


  uint32_t BamRead::sequence_length(void) const
  {
    assert(bam_);
    return bam_->core.l_qseq;
  }


  int32_t BamRead::pos(void) const
  {
    assert(bam_);
    return bam_->core.pos;
  }


  int32_t BamRead::tid(void) const
  {
    assert(bam_);
    return bam_->core.tid;
  }


  int32_t BamRead::mtid(void) const
  {
    assert(bam_);
    return bam_->core.mtid;
  }


  int32_t BamRead::mpos(void) const
  {
    assert(bam_);
    return bam_->core.mpos;
  }


  int32_t BamRead::end(void) const
  {
    assert(bam_);
    return bam_calend(&core(), bam1_cigar(bam_));
  }


  void BamRead::swap(BamRead& other)
  {
    std::swap(bam_, other.bam_);
  }


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

  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
  {
    if (lhs.core().l_qname != rhs.core().l_qname)
      return false;

    return std::equal(lhs.name(),
                      lhs.name() + lhs.core().l_qname,
                      rhs.name());
  }


  bool soft_clipped(const BamRead& bam)
  {
    return left_soft_clipped(bam) || right_soft_clipped(bam);
  }


  uint32_t left_soft_clipped(const BamRead& bam)
  {
    for (size_t i=0; i<bam.core().n_cigar; ++i) {
      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
        return bam.cigar_oplen(i);
      if (bam_cigar_type(bam.cigar_op(i))&1)
        return 0;
    }
    return 0;
  }


  uint32_t right_soft_clipped(const BamRead& bam)
  {
    if (bam.core().n_cigar==0)
      return 0;
    // right clip can't be at first (most left cigar element)
    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
        return bam.cigar_oplen(i);
      if (bam_cigar_type(bam.cigar_op(i))&1)
        return 0;
    }
    return 0;
  }


  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
  {
    return lhs.tid() < rhs.tid() ||
      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
  }


  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
  {
    return lhs.tid() < rhs.tid() ||
      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
  }

}}}