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

#ifndef theplu_yat_omic_bam_read
#define theplu_yat_omic_bam_read

// $Id: BamRead.h 3306 2014-08-21 04:37:21Z peter $

/*
  Copyright (C) 2012, 2013, 2014 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 "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 <functional>
#include <string>
#include <vector>

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

       \deprecated Provided for backward compatibility with 0.11
       API. Use cigar(const utility::Aligner::Cigar&) instead.
    */
    void cigar(const std::vector<uint32_t>& c) YAT_DEPRECATE;

    /**
       \brief set CIGAR

       \param cigar new cigar

       \since new in yat 0.12
    */
    void cigar(const utility::Aligner::Cigar& cigar);

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

       \see Preprocessor defines BAM_F*

       \since implemented since yat 0.12
     */
    uint16_t flag(void) const;

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

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

    /**
       \return query name

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

    /**
       \brief modify name

       \since New in yat 0.11
    */
    void name(const std::string& n);

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

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

    /**
       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);

    /**
       \brief set sequence and quality

       \since New in yat 0.11
     */
    void sequence(const std::string& seq, const std::vector<uint8_t>& qual);

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

    /**
       Exchanging this read with \a other.

       \see swap(BamRead&, BamRead&)
     */
    void swap(BamRead& other);

    /**
       \brief chromosome ID
     */
    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;
    uint32_t calend(const bam1_core_t *c, const uint32_t *cigar) const;
  };

  /**
     Swap specialization for BamRead that is faster than generic
     std::swap as it just swap a pair of 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