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

Last change on this file since 3317 was 3317, checked in by Peter, 7 years ago

refs #806 add test for Pileup class

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.4 KB
Line 
1#ifndef theplu_yat_omic_pileup
2#define theplu_yat_omic_pileup
3
4// $Id: Pileup.h 3317 2014-09-19 05:41:50Z 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  /**
36     \since new in yat 0.13
37   */
38  template<typename Iterator>
39  class Pileup
40  {
41  public:
42
43    /**
44       Essentially a BamRead but with additional information regarding
45       base we currently pointing to
46     */
47    class Entry
48    {
49    public:
50      /**
51         Create an Entry with first read \a b
52       */
53      Entry(const yat::omic::BamRead& b);
54
55      /**
56         \return const ref to BamRead
57       */
58      const yat::omic::BamRead& bam(void) const { return bam_; }
59
60      /**
61         \return cigar element at current position
62       */
63      uint32_t cigar(void) const
64      {
65        assert(c_index_ < bam_.core().n_cigar);
66        return bam_.cigar(c_index_);
67      }
68
69      /**
70         \return base quality of current position
71       */
72      unsigned short qual(void) const
73      {
74        assert(qpos_ < bam_.sequence_length());
75        return bam_.qual(qpos_);
76      }
77
78      /**
79         \return base of current position
80       */
81      uint8_t sequence(void) const;
82
83      /**
84         \brief iterate to next position with coverage
85
86         If any read has an insertion here, next position will be a
87         skip_ref locus; otherwise we iterate to next locus on
88         reference that has coverage.
89       */
90      void increment(void);
91
92      /**
93         \return true if beyond last (aligned) base, i.e., we don't
94         count clipped bases
95       */
96      bool end_of_sequence(void) const
97      {
98        return qpos_ >= bam_.sequence_length();
99      }
100
101      /**
102         \return true if this position was deleted in this query
103       */
104      bool is_del(void) const
105      {
106        assert(c_index_ < bam_.core().n_cigar);
107        return bam_cigar_type(cigar())==2;
108      }
109
110      /**
111         \return index of base at this position
112       */
113      size_t qpos(void) const { return qpos_; }
114    private:
115      yat::omic::BamRead bam_;
116      // index of base pointed to
117      size_t qpos_;
118      size_t c_index_;
119      size_t c_index2_;
120      size_t skip_;
121      void step_cigar(void);
122      uint32_t cigar_op(void) const
123      {
124        return bam_cigar_op(cigar());
125      }
126      friend class Pileup;
127    };
128
129    /**
130       \brief Const iterator that can be used to iterate over reads at
131       current position
132     */
133    typedef typename std::deque<Entry>::const_iterator const_iterator;
134
135    /**
136       Create an empty Pileup
137     */
138    Pileup(void) {}
139
140    /**
141       Create a Pileup that will use reads in range [first, last)
142     */
143    Pileup(Iterator first, Iterator last);
144
145    /**
146       first Entry in pileup
147     */
148    const_iterator begin(void) const { return data_.begin(); }
149
150    /**
151       last Entry in Pileup
152     */
153    const_iterator end(void) const;
154
155    /**
156       \return true if there are more positions to increment into
157     */
158    bool good(void) const { return !(bam_ == bam_end_ && data_.empty()); }
159
160    /**
161       \brief step to next position with coverage
162     */
163    void increment(void);
164
165    /**
166       \return template id as defined bam file header
167     */
168    int32_t tid(void) const { return tid_; }
169
170    /**
171       \return position
172     */
173    int32_t pos(void) const { return pos_; }
174
175    /**
176       \return true if position does not exist in reference, i.e.,
177       inserted region.
178     */
179    bool skip_ref(void) const;
180  private:
181    void step_cigar(void);
182
183    Iterator bam_;
184    Iterator bam_end_;
185    int32_t pos_;
186    int32_t tid_;
187    std::deque<Entry> data_;
188    uint32_t skip_ref_;
189  };
190
191
192  /**
193     \relates Pileup
194
195     \since new in yat 0.13
196   */
197  template<typename Iterator>
198  Pileup<Iterator> make_pileup(Iterator first, Iterator last)
199  {
200    return Pileup<Iterator>(first, last);
201  }
202
203  // implementations
204
205  template<typename Iterator>
206  Pileup<Iterator>::Pileup(Iterator first, Iterator last)
207    : bam_(first), bam_end_(last), skip_ref_(0)
208  {
209    tid_ = 0;
210    pos_ = 0;
211    increment();
212  }
213
214
215  template<typename Iterator>
216  typename Pileup<Iterator>::const_iterator Pileup<Iterator>::end(void) const
217  {
218    if (data_.empty())
219      return data_.end();
220    // if last entry is passed pos, then don't include it in [begin, end) range
221    if ((data_.back().bam().tid() > tid() ||
222         data_.back().bam().pos() > pos())) {
223      return data_.begin() + (data_.size()-1);
224    }
225
226    return data_.end();
227  }
228
229
230  template<typename Iterator>
231  void Pileup<Iterator>::increment(void)
232  {
233    if (data_.empty()) {
234      if (bam_==bam_end_)
235        return;
236      tid_ = bam_->core().tid;
237      pos_ = bam_->core().pos;
238    }
239    else {
240
241      // in insertion
242      if (skip_ref()) {
243        for (size_t i=0; i<data_.size(); ++i) {
244          // ignore reads left of current position
245          if (data_[i].end_of_sequence())
246            continue;
247          if (bam_cigar_type(data_[i].cigar())==1) {
248            data_[i].increment();
249          }
250          --data_[i].skip_;
251        }
252        --skip_ref_;
253      }
254      // walk one base forward on reference
255      else {
256        ++pos_;
257        for (size_t i=0; i<data_.size(); ++i) {
258          if (data_[i].end_of_sequence())
259            continue;
260          data_[i].increment();
261          // check if we're stepping into insertion
262          if (!data_[i].end_of_sequence() && (data_[i].cigar_op() == BAM_CINS))
263            skip_ref_ = std::max(skip_ref_, bam_cigar_oplen(data_[i].cigar()));
264        }
265        // stepped into an insertion
266        if (skip_ref_) {
267          for (size_t i=0; i<data_.size(); ++i)
268            data_[i].skip_ = skip_ref_;
269        }
270      }
271    }
272
273    // read data
274    while (bam_!=bam_end_ && bam_->core().tid==tid_ && bam_->pos()<=pos_) {
275      // skip unmapped reads
276      if (!(bam_->core().flag & BAM_FUNMAP)) {
277        data_.push_back(Entry(*bam_));
278      }
279      ++bam_;
280    }
281
282    // pop reads that are end of sequence
283    while (data_.size() && data_[0].end_of_sequence()) {
284      data_.pop_front();
285    }
286  }
287
288
289  template<typename Iterator>
290  bool Pileup<Iterator>::skip_ref(void) const
291  {
292    return skip_ref_;
293  }
294
295
296  template<typename Iterator>
297  Pileup<Iterator>::Entry::Entry(const yat::omic::BamRead& b)
298    : bam_(b), qpos_(0), c_index_(0), c_index2_(0), skip_(0)
299  {
300    if (bam_.core().n_cigar==0)
301      return;
302    // iterate cigar to first match
303    while (bam_cigar_type(bam_.cigar(c_index_)) != 3) {
304      assert(bam_cigar_type(bam_.cigar(c_index_))!=2);
305      if (bam_cigar_type(bam_.cigar(c_index_)) == 1) {
306        qpos_ += b.cigar_oplen(c_index_);
307      }
308      ++c_index_;
309      if (c_index_ >= bam_.core().n_cigar) {
310        qpos_ = bam_.sequence_length();
311        return;
312      }
313      assert(c_index_ < bam_.core().n_cigar);
314    }
315  }
316
317
318  template<typename Iterator>
319  void Pileup<Iterator>::Entry::increment(void)
320  {
321    if (bam_cigar_type(cigar()) & 1)
322      ++qpos_;
323    step_cigar();
324  }
325
326
327  template<typename Iterator>
328  uint8_t Pileup<Iterator>::Entry::sequence(void) const
329  {
330    if (skip_==0) {
331      if (bam_cigar_type(cigar()) & 1)
332        return bam_.sequence(qpos_);
333      else
334        return 0;
335    }
336    if (cigar_op()!=BAM_CINS)
337      return 0;
338    return bam_.sequence(qpos_);
339  }
340
341
342  template<typename Iterator>
343  void Pileup<Iterator>::Entry::step_cigar(void)
344  {
345    ++c_index2_;
346    assert(c_index_<bam_.core().n_cigar);
347    if (c_index2_ >= bam_.cigar_oplen(c_index_)) {
348      c_index2_ = 0;
349      ++c_index_;
350    }
351  }
352
353}}}
354#endif
Note: See TracBrowser for help on using the repository browser.