Pileup.h

#ifndef theplu_yat_omic_pileup
#define theplu_yat_omic_pileup

// $Id: Pileup.h 3855 2020-01-02 01:11:34Z peter $

/*
  Copyright (C) 2014, 2016, 2017, 2018, 2019 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "BamRead.h"

#include "yat/utility/CigarIterator.h"
#include "yat/utility/yat_assert.h"

#include <boost/iterator/iterator_traits.hpp>
#include <boost/iterator/iterator_concepts.hpp>
#include <boost/concept/assert.hpp>
#include <boost/shared_ptr.hpp>

#include <algorithm>
#include <list>

namespace theplu {
namespace yat {
namespace omic {

  template<typename Iterator>
  class Pileup
  {
  public:
    class Entry
    {
    public:
      Entry(void) {}

      Entry(const BamRead& b);

      const BamRead& bam(void) const
      {
        YAT_ASSERT(bam_.get());
        return *bam_;
      }

      uint8_t cigar_op(void) const { return *cigar_; }


      uint32_t cigar_oplen(void) const
      {
        YAT_ASSERT(bam_.get());
        YAT_ASSERT(cigar_.base() >= bam_->cigar());
        return bam_cigar_oplen(*cigar_.base());
      }


      unsigned short qual(void) const
      {
        YAT_ASSERT(bam_.get());
        YAT_ASSERT(qpos_ < bam_->sequence_length());
        return bam_->qual(qpos_);
      }

      uint8_t sequence(void) const;

      void increment(void);

      bool is_del(void) const
      {
        return cigar_type()==2;
      }

      size_t qpos(void) const { return qpos_; }
    private:
      boost::shared_ptr<BamRead> bam_;
      // index of base pointed to
      size_t qpos_;
      yat::utility::CigarIterator<const uint32_t*> cigar_;
      size_t skip_;
      // return 0, 1, 2, or 3 see bam_cigar_type in yat/utility/Cigar.h
      uint8_t cigar_type(void) const { return bam_cigar_type(cigar_op()); }
      friend class Pileup;
    };

  public:
    typedef typename std::list<Entry>::const_iterator const_iterator;

    typedef typename std::list<Entry>::const_reverse_iterator
    const_reverse_iterator;

    Pileup(void) {}

    Pileup(Iterator first, Iterator last);

    const_iterator begin(void) const;

    const_iterator end(void) const;

    bool good(void) const { return !(bam_ == bam_end_ && data_.empty()); }

    void increment(void);

    size_t size(void) const;

    const_reverse_iterator rbegin(void) const;

    const_reverse_iterator rend(void) const;

    int32_t tid(void) const { return tid_; }

    int32_t pos(void) const { return pos_; }

    bool skip_ref(void) const;
  private:
    typedef typename std::list<Entry>::iterator iterator;

    Iterator bam_;
    Iterator bam_end_;
    int32_t pos_;
    int32_t tid_;
    std::list<Entry> data_;
    uint32_t skip_ref_;
  };


  template<typename Iterator>
  Pileup<Iterator> make_pileup(Iterator first, Iterator last)
  {
    return Pileup<Iterator>(first, last);
  }

  // implementations

  template<typename Iterator>
  Pileup<Iterator>::Pileup(Iterator first, Iterator last)
    : bam_(first), bam_end_(last), pos_(0), tid_(0), skip_ref_(0)
  {
    BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<Iterator>));
    BOOST_CONCEPT_ASSERT((boost_concepts::SinglePassIterator<Iterator>));
    using boost::Convertible;
    typedef typename boost::iterator_value<Iterator>::type value_type;
    BOOST_CONCEPT_ASSERT((Convertible<value_type, BamRead>));
    increment();
  }


  template<typename Iterator>
  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::begin(void) const
  {
    return data_.begin();
  }


  template<typename Iterator>
  typename Pileup<Iterator>::const_reverse_iterator
  Pileup<Iterator>::rbegin(void) const
  {
    return data_.rbegin();
  }


  template<typename Iterator>
  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
  {
    return data_.end();
  }


  template<typename Iterator>
  typename Pileup<Iterator>::const_reverse_iterator
  Pileup<Iterator>::rend(void) const
  {
    return data_.rend();
  }


  template<typename Iterator>
  void Pileup<Iterator>::increment(void)
  {
    if (!data_.empty()) {
      // in insertion
      if (skip_ref()) {
        for (iterator it = data_.begin(); it!=data_.end(); ++it) {
          if (it->cigar_type() == 1) {
            it->increment();
          }
          --it->skip_;
        }
        --skip_ref_;
      }
      // walk one base forward on reference
      else {
        ++pos_;
        typename std::list<Entry>::iterator read = data_.begin();
        while (read != data_.end()) {
          read->increment();
          // remove entry if read is entirely left of pos
          if (read->bam().end() <= pos_) { // FIXME speedup possible
            data_.erase(read++);
            continue;
          }
          // check if we're stepping into insertion
          if (read->cigar_op() == BAM_CINS) {
            skip_ref_ = std::max(skip_ref_, read->cigar_oplen());
          }
          ++read;
        }

        // stepped into an insertion
        if (skip_ref_)
          for (iterator iter=data_.begin(); iter!=data_.end(); ++iter)
            iter->skip_ = skip_ref_;
      }
    }


    // fast forward to next covered locus
    if (data_.empty()) {
      if (bam_==bam_end_)
        return;
      tid_ = bam_->core().tid;
      pos_ = bam_->core().pos;
    }

    // read data
    while (bam_!=bam_end_ && bam_->core().tid==tid_ && bam_->pos()<=pos_) {
      // skip unmapped reads
      if (!(bam_->core().flag & BAM_FUNMAP)) {
        data_.push_back(Entry(*bam_));
      }
      ++bam_;
    }
  }


  template<typename Iterator>
  bool Pileup<Iterator>::skip_ref(void) const
  {
    return skip_ref_;
  }


  template<typename Iterator>
  size_t Pileup<Iterator>::size(void) const
  {
    return data_.size();
  }


  template<typename Iterator>
  Pileup<Iterator>::Entry::Entry(const BamRead& b)
    : bam_(new BamRead(b)), qpos_(0), cigar_(bam_->cigar()), skip_(0)
  {
    YAT_ASSERT(cigar_.base() == bam_->cigar());
    if (bam_->core().n_cigar==0)
      return;
    // iterate cigar to first match
    while (qpos_ < bam_->sequence_length() && cigar_type() != 3) {
      // type is either 0 (hardclipped) or 1 (softclipped)
      YAT_ASSERT(cigar_type()!=2);
      if (cigar_type() == 1)
        ++qpos_;
      ++cigar_;
    }
  }


  template<typename Iterator>
  void Pileup<Iterator>::Entry::increment(void)
  {
    YAT_ASSERT(cigar_.base() >= bam_->cigar());
    if (cigar_type() & 1)
      ++qpos_;
    ++cigar_;
    YAT_ASSERT(cigar_.base() >= bam_->cigar());
  }


  template<typename Iterator>
  uint8_t Pileup<Iterator>::Entry::sequence(void) const
  {
    if (skip_==0) {
      if (cigar_type() & 1)
        return bam_->sequence(qpos_);
      return 0;
    }
    if (cigar_op()!=BAM_CINS)
      return 0;
    return bam_->sequence(qpos_);
  }

}}}
#endif