source: trunk/yat/omic/BamRead.h @ 2985

#ifndef theplu_yat_omic_bam_read
#define theplu_yat_omic_bam_read

// $Id: BamRead.h 2985 2013-02-18 07:36:14Z peter $
//
// Copyright (C) 2012, 2013 Peter Johansson
//
// This program 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.
//
// This program 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_BAM_HEADER
#include YAT_SAM_HEADER

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

// backport #defines from samtools 0.1.19
#ifndef BAM_CIGAR_STR
#define BAM_CIGAR_STR "MIDNSHP=XB"
// lookup table used in bam_cigar_type
#define BAM_CIGAR_TYPE 0x3C1A7

/// index of operation in range [0, 9], see below.
#define bam_cigar_op(c) ((c)&BAM_CIGAR_MASK)
// length of the operation
#define bam_cigar_oplen(c) ((c)>>BAM_CIGAR_SHIFT)
// char describing the operation, see column 2 in table below
#define bam_cigar_opchr(c) (BAM_CIGAR_STR[bam_cigar_op(c)])
// Create a cigar element with length l and operation o
#define bam_cigar_gen(l, o) ((l)<<BAM_CIGAR_SHIFT|(o))
// Returns a two-bits number describing the type of operation. The
// first bit, bam_cigar_type & 1, (4th column in table below) is 1 iff
// operation consumes query sequence. The second bit, bam_cigar_type &
// 2, (3rd column below) is 1 iff operation consumes reference
// sequence.
#define bam_cigar_type(o) (BAM_CIGAR_TYPE>>((o)<<1)&3)
#endif

/*
  CIGAR
0 M 1 1 match
1 I 0 1 insert
2 D 1 0 deletion
3 N 1 0 spliced (rna)
4 S 0 1 soft clipped
5 H 0 0 hard clipped
6 P 0 0 padded
7 = 1 1 equal
8 X 1 1 mismatch
9 B 0 0

3rd column: consumes reference sequence
4th column: consumes query sequence
 */

namespace theplu {
namespace yat {
namespace omic {

  /**
     \brief Class holding a bam query

     This is a wrapper around bam1_t and most of its information is
     held in the core struct. A BamRead is typically created from a
     InBamFile.

     \see samtools

     \since New in yat 0.10
   */
  class BamRead
  {
  public:
    /**
       \brief default constructor

       Constructed object contains no data and most operations are not
       defined
     */
    BamRead(void);

    /**
       \brief Copy constructor
     */
    BamRead(const BamRead& other);

    /**
       \brief Destructor
     */
    virtual ~BamRead(void);

    /**
       \brief assignment operator
     */
    const BamRead& operator=(const BamRead& rhs);

    /**
       \return pointer to auxiliary data
     */
    const uint8_t* aux(void) const;

    /**
       \brief access core data struct

       \see samtools C api documentaion
     */
    const bam1_core_t& core(void) const;

    /**
       \brief access core data struct

       \see samtools C api documentaion
     */
    bam1_core_t& core(void);

    /**
       \brief access CIGAR array

       In each element the lower 4 bits gives a CIGAR operation and
       the upper 28 bits keep the length.
     */
    const uint32_t* cigar(void) const;

    /**
       \return \a i th element of CIGAR array
     */
    uint32_t cigar(size_t i) const;

    /**
       \return \a i th CIGAR operation
     */
    uint32_t cigar_op(size_t i) const;

    /**
       \return length of \a i th CIGAR element
     */
    uint32_t cigar_oplen(size_t i) const;

    /**
       Translate CIGAR array to a string such as '72M3S'
     */
    std::string cigar_str(void) const;

    /**
       \brief set CIGAR

       \param c new cigar
    */
    void cigar(const std::vector<uint32_t>& c);

    /**
       \return Quality of base \a i
     */
    uint8_t qual(size_t i) const;

    /**
       \return query name

       Length of array is described by core().l_qname
     */
    const char* name(void) const;

    /**
       \brief chromosome ID
     */
    int32_t tid(void) const;

    /**
       \brief 0-based laftmost coordinate
     */
    int32_t pos(void) const;

    /**
       \brief rightmost coordinate

       Coordinate is 0-based, i.e., the end is one passed the last
       matching position.

       \see bam_calend
     */
    int32_t end(void) const;

    /**
       \brief Chromosome ID for mate
     */
    int32_t mtid(void) const;

    /**
       \brief leftmost position for mate
     */
    int32_t mpos(void) const;

    /**
       Each character in returned sequence is either A, C, G, T, or N.

       \return sequence
     */
    std::string sequence(void) const;

    /**
       4-bit integer describing base \a index

       \see bam_nt16_rev_table
     */
    uint8_t sequence(size_t index) const;

    /**
       \brief modify a base in sequence

       Set i-th base in sequence to \a x, where seq is a 4-bit integer.

       \see bam_nt16_table
     */
    void sequence(size_t i, uint8_t x);

    /**
       \see core().l_qseq
     */
    uint32_t sequence_length(void) const;

    /**
       \brief bitwise flag

       \see Preprocessor defines BAM_F*
     */
    uint16_t flag(void) const;

    /**
       Exchanging this read with \a other.
     */
    void swap(BamRead& other);

  private:
    bam1_t* bam_;

    friend class InBamFile;
    friend class OutBamFile;
    friend class BamReadIterator;
  };

  /**
     Swap specialization for BamRead that is faster than generic
     std::swap as it just swap a pointers.

     \since New in yat 0.10
   */
  void swap(BamRead& lhs, BamRead& rhs);

  /**
     \return \c true if read is soft clipped, either left_soft_clipped
     or right_soft_clipped.

     \since New in yat 0.10
   */
  bool soft_clipped(const BamRead& bam);

  /**
     If read is soft clipped on left side, return how many bases are
     clipped, otherwise return 0.

     \since New in yat 0.10
   */
  uint32_t left_soft_clipped(const BamRead& bam);

  /**
     If read is soft clipped on right side, return how many bases are
     clipped, otherwise return 0.

     \since New in yat 0.10
   */
  uint32_t right_soft_clipped(const BamRead& bam);

  /**
     return \c true if query names are equal

     \see BamRead::name()

     \since New in yat 0.10
   */
  bool same_query_name(const BamRead& lhs, const BamRead& rhs);

  /**
     Functor to compare two reads with respect to their leftmost
     coordinate.

     \see BamRead

     \since New in yat 0.10
   */
  struct BamLessPos
    : public std::binary_function<const BamRead&, const BamRead&, bool>
  {
    /**
       \return true if lhs tid is less than rhs tid; or tid
       are equal and lhs pos is smaller than rhs pos.

       \see BamRead::tid() and BamRead::pos()
     */
    bool operator()(const BamRead& lhs, const BamRead& rhs) const;
  };


  /**
     Functor to compare two reads with respect to their rightmost
     coordinate.

     \see BamRead

     \since New in yat 0.10
   */
  struct BamLessEnd
    : public std::binary_function<const BamRead&, const BamRead&, bool>
  {
    /**
       \return true if lhs tid is less than rhs tid; or tid
       are equal and lhs end is smaller than rhs end.

       \see BamRead::tid() and BamRead::end()
     */
    bool operator()(const BamRead& lhs, const BamRead& rhs) const;
  };

}}}
#endif