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

#ifndef theplu_yat_omic_bam_read
#define theplu_yat_omic_bam_read

// $Id: BamRead.h 3029 2013-04-21 01:04:41Z 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_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

       \see aux_size(void)
     */
    const uint8_t* aux(void) const;

    /**
       Use bam_aux2? functions (in samtools C api) to convert returned
       pointer to corresponding type.

       \return pointer to field associated with \a tag, NULL if \a tag
       doesn't exist.

       \see bam_aux_get

       \since New in yat 0.11
     */
    const uint8_t* aux(const char tag[2]) const;

    /**
       \brief append a new tag to aux field

       \param tag two-charcter tag to append
       \param type describes which type and can be 'iIsScCdfAZH'
       \param len length of data
       \param data pointer to data

       \since New in yat 0.11

       \see SAM specification
     */
    void aux_append(const char tag[2], char type, int len, uint8_t* data);

    /**
       \brief remove a tag in aux field

       \throw utility::runtime_error if \a tag is not present in read

       \since New in yat 0.11
     */
    void aux_del(const char tag[2]);

    /**
       \return length of aux field

       \since New in yat 0.11
     */
    int aux_size(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;

    /**
       \brief set quality of a base

       \param i base to modify
       \param q new quality

       \since New in yat 0.11
     */
    void qual(size_t i, uint8_t q);

    /**
       \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

       \since New in yat 0.11
     */
    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

     \relates BamRead
   */
  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

     \relates BamRead
   */
  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

     \relates BamRead
   */
  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

     \relates BamRead
   */
  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

     \relates BamRead
   */
  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