BamRead.h

#ifndef theplu_yat_omic_bam_read
#define theplu_yat_omic_bam_read

// $Id: BamRead.h 3550 2017-01-03 05:41:02Z peter $

/*
  Copyright (C) 2012, 2013, 2014, 2015, 2016 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_bam.h"

#include "yat/utility/Aligner.h"
// This file has to be included to keep compatibility with yat 0.11
#include "yat/utility/Cigar.h"
#include "yat/utility/deprecate.h"

#include YAT_SAM_HEADER

#include <functional>
#include <string>
#include <vector>

namespace theplu {
namespace yat {
namespace omic {

  class BamRead
  {
  public:
    BamRead(void);

    BamRead(const BamRead& other);

    virtual ~BamRead(void);

    const BamRead& operator=(const BamRead& rhs);

    const uint8_t* aux(void) const;

    const uint8_t* aux(const char tag[2]) const;

    void aux_append(const char tag[2], char type, int len, uint8_t* data);

    void aux_del(const char tag[2]);

    int aux_size(void) const;

    const bam1_core_t& core(void) const;

    bam1_core_t& core(void);

    const uint32_t* cigar(void) const;

    uint32_t cigar(size_t i) const;

    uint32_t cigar_op(size_t i) const;

    uint32_t cigar_oplen(size_t i) const;

    std::string cigar_str(void) const;

    void cigar(const std::vector<uint32_t>& c) YAT_DEPRECATE;

    void cigar(const utility::Aligner::Cigar& cigar);

    int32_t end(void) const;

    uint16_t flag(void) const;

    int32_t mpos(void) const;

    int32_t mtid(void) const;

    const char* name(void) const;

    void name(const std::string& n);

    int32_t pos(void) const;

    uint8_t qual(size_t i) const;

    void qual(size_t i, uint8_t q);

    std::string sequence(void) const;

    uint8_t sequence(size_t index) const;

    char sequence_str(size_t index) const;

    void sequence(size_t i, uint8_t x);

    void sequence(const std::string& seq, const std::vector<uint8_t>& qual);

    uint32_t sequence_length(void) const;

    void swap(BamRead& other);

    int32_t tid(void) const;

  private:
    // ensure capacity of data pointer is (at least) n
    void reserve(int n);

    bam1_t* bam_;

    friend class InBamFile;
    friend class OutBamFile;
    friend class BamReadIterator;
#ifndef YAT_HAVE_HTSLIB
    uint32_t calend(const bam1_core_t *c, const uint32_t *cigar) const;
#endif
    // access data length, current length of data
    int& data_size(void);
    const int& data_size(void) const;
  };

  const unsigned char nt16_table(char);

  const char nt16_str(uint8_t);

  void swap(BamRead& lhs, BamRead& rhs);

  bool soft_clipped(const BamRead& bam);

  uint32_t left_soft_clipped(const BamRead& bam);

  uint32_t right_soft_clipped(const BamRead& bam);

  bool same_query_name(const BamRead& lhs, const BamRead& rhs);

  bool less_query_name(const BamRead& lhs, const BamRead& rhs);

  struct BamLessPos
    : public std::binary_function<const BamRead&, const BamRead&, bool>
  {
    bool operator()(const BamRead& lhs, const BamRead& rhs) const;
  };


  struct BamLessEnd
    : public std::binary_function<const BamRead&, const BamRead&, bool>
  {
    bool operator()(const BamRead& lhs, const BamRead& rhs) const;
  };

  struct BamLessName
    : public std::binary_function<const BamRead&, const BamRead&, bool>
  {
    bool operator()(const BamRead& lhs, const BamRead& rhs) const;
  };


  template<typename RandomAccessIterator>
  double alignment_score(const BamRead& bam, RandomAccessIterator ref,
                         double mismatch, double gap=0.0, double open_gap=0.0)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<RandomAccessIterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<RandomAccessIterator>));
    double res = 0;
    size_t q = 0;

    for (uint32_t i = 0; i<bam.core().n_cigar; ++i) {
      switch (bam.cigar_op(i)) {
      case (BAM_CMATCH): {
        RandomAccessIterator end = ref + bam.cigar_oplen(i);
        for (; ref!=end; ++ref) {
          if (bam.sequence_str(q) == *ref)
            res += 1.0;
          else
            res -= mismatch;
          ++q;
        }
        break;
      }
      case (BAM_CEQUAL): {
        q += bam.cigar_oplen(i);
        ref += bam.cigar_oplen(i);
        res += bam.cigar_oplen(i);
        break;
      }
      case (BAM_CDIFF): {
        q += bam.cigar_oplen(i);
        ref += bam.cigar_oplen(i);
        res -= bam.cigar_oplen(i)*mismatch;
        break;
      }
      case (BAM_CSOFT_CLIP): {
        q += bam.cigar_oplen(i);
        break;
      }
      case (BAM_CINS): {
        q += bam.cigar_oplen(i);
        res -= open_gap + bam.cigar_oplen(i) * gap;
        break;
      }
      case (BAM_CDEL): {
        ref += bam.cigar_oplen(i);
        res -= open_gap + bam.cigar_oplen(i) * gap;
        break;
      }
      } // end switch
    }
    return res;
  }

}}}
#endif