source: trunk/yat/omic/Pileup.h @ 3310

Last change on this file since 3310 was 3310, checked in by Peter, 8 years ago

a first version of new class Pileup. refs #806

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.2 KB
Line 
1#ifndef theplu_yat_omic_pileup
2#define theplu_yat_omic_pileup
3
4// $Id: Pileup.h 3310 2014-08-25 06:52:19Z peter $
5
6/*
7  Copyright (C) 2012, 2013, 2014 Peter Johansson
8
9  This file is part of the yat library, http://dev.thep.lu.se/yat
10
11  The yat library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU General Public License as
13  published by the Free Software Foundation; either version 3 of the
14  License, or (at your option) any later version.
15
16  The yat library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program. If not, see <http://www.gnu.org/licenses/>.
23*/
24
25#include "BamRead.h"
26
27#include <algorithm>
28#include <cassert>
29#include <deque>
30
31namespace theplu {
32namespace yat {
33namespace omic {
34
35  template<typename Iterator>
36  class Pileup
37  {
38  public:
39
40    /**
41       Essentially a BamRead but with additional information regarding
42       base we currently pointing to
43     */
44    class Entry
45    {
46    public:
47      /**
48         Create an Entry with first read \a b
49       */
50      Entry(const yat::omic::BamRead& b);
51
52      /**
53         \return const ref to BamRead
54       */
55      const yat::omic::BamRead& bam(void) const { return bam_; }
56
57      uint32_t cigar(void) const
58      {
59        assert(c_index_ < bam_.core().n_cigar);
60        return bam_.cigar(c_index_);
61      }
62
63      /**
64         \return base quality of current position
65       */
66      unsigned short qual(void) const
67      {
68        assert(qpos_ < bam_.sequence_length());
69        return bam_.qual(qpos_);
70      }
71
72      /**
73         \return base of current position
74       */
75      uint8_t sequence(void) const;
76
77      /**
78         \brief iterate to next position with coverage
79
80         If any read has an insertion here, next position will be a
81         skip_ref locus; otherwise we iterate to next locus on
82         reference that has coverage.
83       */
84      void increment(void);
85
86      /**
87         \return true if beyond last (aligned) base, i.e., we don't
88         count clipped bases
89       */
90      bool end_of_sequence(void) const
91      {
92        return qpos_ >= bam_.sequence_length();
93      }
94
95      /**
96         \return true if this position was deleted in this query
97       */
98      bool is_del(void) const
99      {
100        assert(c_index_ < bam_.core().n_cigar);
101        return bam_cigar_type(cigar())==2;
102      }
103
104      size_t qpos(void) const { return qpos_; }
105    private:
106      yat::omic::BamRead bam_;
107      // index of base pointed to
108      size_t qpos_;
109      size_t c_index_;
110      size_t c_index2_;
111      size_t skip_;
112      void step_cigar(void);
113      uint32_t cigar_op(void) const
114      {
115        return bam_cigar_op(cigar());
116      }
117      friend class Pileup;
118    };
119
120    // FIXME change to a filter iterator that ignores end_of_sequence entries
121    typedef typename std::deque<Entry>::const_iterator const_iterator;
122
123    /**
124       Create an empty Pileup
125     */
126    Pileup(void) {}
127
128    /**
129       Create a Pileup that will use reads in range [first, last)
130     */
131    Pileup(Iterator first, Iterator last);
132
133    /**
134       first Entry in pileup
135     */
136    const_iterator begin(void) const { return data_.begin(); }
137
138    /**
139       last Entry in Pileup
140     */
141    const_iterator end(void) const;
142
143    /**
144       \return true if there are more positions to increment into
145     */
146    bool good(void) const { return !(bam_ == bam_end_ && data_.empty()); }
147
148    /**
149     */
150    void increment(void);
151
152    /**
153       \return template id as defined bam file header
154     */
155    int32_t tid(void) const { return tid_; }
156
157    /**
158       \return position
159     */
160    int32_t pos(void) const { return pos_; }
161
162    /**
163       \return true if position does not exist in reference, i.e.,
164       inserted region.
165     */
166    bool skip_ref(void) const;
167  private:
168    void step_cigar(void);
169
170    Iterator bam_;
171    Iterator bam_end_;
172    int32_t pos_;
173    int32_t tid_;
174    std::deque<Entry> data_;
175    uint32_t skip_ref_;
176  };
177
178
179  template<typename Iterator>
180  Pileup<Iterator> make_pileup(Iterator first, Iterator last)
181  {
182    return Pileup<Iterator>(first, last);
183  }
184
185  // implementations
186
187  template<typename Iterator>
188  Pileup<Iterator>::Pileup(Iterator first, Iterator last)
189    : bam_(first), bam_end_(last), skip_ref_(0)
190  {
191    tid_ = 0;
192    pos_ = 0;
193    increment();
194  }
195
196
197  template<typename Iterator>
198  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
199  {
200    if (data_.empty())
201      return data_.end();
202    // if last entry is passed pos, then don't include it in [begin, end) range
203    if ((data_.back().bam().tid() > tid() ||
204         data_.back().bam().pos() > pos())) {
205      return data_.begin() + (data_.size()-1);
206    }
207
208    return data_.end();
209  }
210
211
212  template<typename Iterator>
213  void Pileup<Iterator>::increment(void)
214  {
215    if (data_.empty()) {
216      if (bam_==bam_end_)
217        return;
218      tid_ = bam_->core().tid;
219      pos_ = bam_->core().pos;
220    }
221    else {
222
223      // in insertion
224      if (skip_ref()) {
225        for (size_t i=0; i<data_.size(); ++i) {
226          // ignore reads left of current position
227          if (data_[i].end_of_sequence())
228            continue;
229          if (bam_cigar_type(data_[i].cigar())==1) {
230            data_[i].increment();
231          }
232          --data_[i].skip_;
233        }
234        --skip_ref_;
235      }
236      // walk one base forward on reference
237      else {
238        ++pos_;
239        for (size_t i=0; i<data_.size(); ++i) {
240          if (data_[i].end_of_sequence())
241            continue;
242          data_[i].increment();
243          // check if we're stepping into insertion
244          if (!data_[i].end_of_sequence() && (data_[i].cigar_op() == BAM_CINS))
245            skip_ref_ = std::max(skip_ref_, bam_cigar_oplen(data_[i].cigar()));
246        }
247        // stepped into an insertion
248        if (skip_ref_) {
249          for (size_t i=0; i<data_.size(); ++i)
250            data_[i].skip_ = skip_ref_;
251        }
252      }
253    }
254
255    // read data
256    while (bam_!=bam_end_ && bam_->core().tid==tid_ && bam_->pos()<=pos_) {
257      // skip unmapped reads
258      if (!(bam_->core().flag & BAM_FUNMAP)) {
259        data_.push_back(Entry(*bam_));
260      }
261      ++bam_;
262    }
263
264    // pop reads that are end of sequence
265    while (data_.size() && data_[0].end_of_sequence()) {
266      data_.pop_front();
267    }
268  }
269
270
271  template<typename Iterator>
272  bool Pileup<Iterator>::skip_ref(void) const
273  {
274    return skip_ref_;
275  }
276
277
278  template<typename Iterator>
279  Pileup<Iterator>::Entry::Entry(const yat::omic::BamRead& b)
280    : bam_(b), qpos_(0), c_index_(0), c_index2_(0), skip_(0)
281  {
282    if (bam_.core().n_cigar==0)
283      return;
284    // iterate cigar to first match
285    while (bam_cigar_type(bam_.cigar(c_index_)) != 3) {
286      assert(bam_cigar_type(bam_.cigar(c_index_))!=2);
287      if (bam_cigar_type(bam_.cigar(c_index_)) == 1) {
288        qpos_ += b.cigar_oplen(c_index_);
289      }
290      ++c_index_;
291      if (c_index_ >= bam_.core().n_cigar) {
292        qpos_ = bam_.sequence_length();
293        return;
294      }
295      assert(c_index_ < bam_.core().n_cigar);
296    }
297  }
298
299
300  template<typename Iterator>
301  void Pileup<Iterator>::Entry::increment(void)
302  {
303    if (bam_cigar_type(cigar()) & 1)
304      ++qpos_;
305    step_cigar();
306  }
307
308
309  template<typename Iterator>
310  uint8_t Pileup<Iterator>::Entry::sequence(void) const
311  {
312    if (skip_==0) {
313      if (bam_cigar_type(cigar()) & 1)
314        return bam_.sequence(qpos_);
315      else
316        return 0;
317    }
318    if (cigar_op()!=BAM_CINS)
319      return 0;
320    return bam_.sequence(qpos_);
321  }
322
323
324  template<typename Iterator>
325  void Pileup<Iterator>::Entry::step_cigar(void)
326  {
327    ++c_index2_;
328    assert(c_index_<bam_.core().n_cigar);
329    if (c_index2_ >= bam_.cigar_oplen(c_index_)) {
330      c_index2_ = 0;
331      ++c_index_;
332    }
333  }
334
335}}}
336#endif
Note: See TracBrowser for help on using the repository browser.