source: branches/0.10-stable/yat/omic/BamRead.cc @ 2928

Last change on this file since 2928 was 2928, checked in by Peter, 9 years ago

Support different styles of including bam header files: <bam.h>,
<bam/bam.h>, and <samtools/bam.h>. Style is identified at configure
time and stored in generated 'config_public.h' by defining either
HAVE_BAM_H, HAVE_BAM_BAM_H or HAVE_SAMTOOLS_BAM_H. For convenience,
macros YAT_[B|S]AM_HEADER are defined in 'config_bam.h' that can be
used as '#include YAT_[B|S]AB_HEADER' rather than '#include
<[b|s]am.h'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 5.2 KB
Line 
1// $Id: BamRead.cc 2928 2012-12-25 23:20:10Z peter $
2//
3// Copyright (C) 2012 Peter Johansson
4//
5// This program is free software; you can redistribute it and/or modify
6// it under the terms of the GNU General Public License as published by
7// the Free Software Foundation; either version 3 of the License, or
8// (at your option) any later version.
9//
10// This program is distributed in the hope that it will be useful, but
11// WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13// General Public License for more details.
14//
15// You should have received a copy of the GNU General Public License
16// along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18#include <config.h>
19
20#include "BamRead.h"
21#include "config_bam.h"
22
23#include YAT_BAM_HEADER
24
25#include <algorithm>
26#include <cassert>
27#include <sstream>
28
29namespace theplu {
30namespace yat {
31namespace omic {
32
33  BamRead::BamRead(void)
34    : bam_(bam_init1())
35  {}
36
37
38  BamRead::BamRead(const BamRead& other)
39    : bam_(bam_dup1(other.bam_))
40  {
41  }
42
43
44  BamRead::~BamRead(void)
45  {
46    bam_destroy1(bam_);
47  }
48
49
50  const BamRead& BamRead::operator=(const BamRead& rhs)
51  {
52    if (bam_!=rhs.bam_)
53      bam_copy1(bam_, rhs.bam_);
54    return *this;
55  }
56
57
58  const uint8_t* BamRead::aux(void) const
59  {
60    return bam1_aux(bam_);
61  }
62
63
64  const uint32_t* BamRead::cigar(void) const
65  {
66    return bam1_cigar(bam_);
67  }
68
69
70  uint32_t BamRead::cigar(size_t i) const
71  {
72    assert(i<core().n_cigar);
73    return bam1_cigar(bam_)[i];
74  }
75
76
77  uint32_t BamRead::cigar_op(size_t i) const
78  {
79    assert(i<core().n_cigar);
80    return bam_cigar_op(cigar(i));
81  }
82
83
84  uint32_t BamRead::cigar_oplen(size_t i) const
85  {
86    assert(i<core().n_cigar);
87    return bam_cigar_oplen(cigar(i));
88  }
89
90
91  std::string BamRead::cigar_str(void) const
92  {
93    std::ostringstream os;
94    std::string str = BAM_CIGAR_STR;
95    for (uint32_t i=0; i<core().n_cigar; ++i) {
96      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
97         << str[bam_cigar_op(bam1_cigar(bam_)[i])];
98    }
99    return os.str();
100  }
101
102
103  void BamRead::cigar(const std::vector<uint32_t>& c)
104  {
105    int m_data = bam_->m_data;
106    int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
107    // double capacity if needed
108    if (m_data < data_len) {
109      m_data = data_len;
110      kroundup32(m_data);
111    }
112    uint8_t* data = (uint8_t*)calloc(m_data, 1);
113    /*
114      data comes as qname, cigar, rest
115     */
116    // first copy qname
117    memcpy(data, bam_->data, core().l_qname);
118    // copy new cigar
119    if (c.size())
120      // each cigar element is 4 bytes
121      memcpy(data + core().l_qname, &c[0], 4*c.size());
122    // copy remaining data
123    memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
124           data_len - 4*c.size() - core().l_qname);
125    bam_->m_data = m_data;
126    bam_->data_len = data_len;
127    core().n_cigar = c.size();
128    free(bam_->data);
129    bam_->data = data;
130  }
131
132
133  const bam1_core_t& BamRead::core(void) const
134  { return bam_->core; }
135
136
137  bam1_core_t& BamRead::core(void)
138  { return bam_->core; }
139
140
141  const char* BamRead::name(void) const
142  { return bam1_qname(bam_); }
143
144
145  uint8_t BamRead::qual(size_t i) const
146  { return *(bam1_qual(bam_)+i); }
147
148
149  std::string BamRead::sequence(void) const
150  {
151    // FIXME bam_nt16_rev_table not available in old lbam
152    std::string result(sequence_length(), ' ');
153    for (size_t i=0; i<result.size(); ++i)
154      result[i] = bam_nt16_rev_table[sequence(i)];
155    return result;
156  }
157
158
159  uint8_t BamRead::sequence(size_t index) const
160  {
161    assert(bam_);
162    return bam1_seqi(bam1_seq(bam_),index);
163  }
164
165
166  uint32_t BamRead::sequence_length(void) const
167  {
168    assert(bam_);
169    return bam_->core.l_qseq;
170  }
171
172
173  int32_t BamRead::pos(void) const
174  {
175    assert(bam_);
176    return bam_->core.pos;
177  }
178
179
180  int32_t BamRead::tid(void) const
181  {
182    assert(bam_);
183    return bam_->core.tid;
184  }
185
186
187
188  int32_t BamRead::mpos(void) const
189  {
190    assert(bam_);
191    return bam_->core.mpos;
192  }
193
194
195  int32_t BamRead::end(void) const
196  {
197    assert(bam_);
198    return bam_calend(&core(), bam1_cigar(bam_));
199  }
200
201
202  void BamRead::swap(BamRead& other)
203  {
204    std::swap(bam_, other.bam_);
205  }
206
207
208  void swap(BamRead& lhs, BamRead& rhs)
209  {
210    lhs.swap(rhs);
211  }
212
213  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
214  {
215    if (lhs.core().l_qname != rhs.core().l_qname)
216      return false;
217
218    return std::equal(lhs.name(),
219                      lhs.name() + lhs.core().l_qname,
220                      rhs.name());
221  }
222
223
224  bool soft_clipped(const BamRead& bam)
225  {
226    return left_soft_clipped(bam) || right_soft_clipped(bam);
227  }
228
229
230  uint32_t left_soft_clipped(const BamRead& bam)
231  {
232    for (size_t i=0; i<bam.core().n_cigar; ++i) {
233      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
234        return bam.cigar_oplen(i);
235      if (bam_cigar_type(bam.cigar_op(i))&1)
236        return 0;
237    }
238    return 0;
239  }
240
241
242  uint32_t right_soft_clipped(const BamRead& bam)
243  {
244    if (bam.core().n_cigar==0)
245      return 0;
246    // right clip can't be at first (most left cigar element)
247    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
248      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
249        return bam.cigar_oplen(i);
250      if (bam_cigar_type(bam.cigar_op(i))&1)
251        return 0;
252    }
253    return 0;
254  }
255
256
257  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
258  {
259    return lhs.tid() < rhs.tid() ||
260      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
261  }
262
263
264  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
265  {
266    return lhs.tid() < rhs.tid() ||
267      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
268  }
269
270}}}
Note: See TracBrowser for help on using the repository browser.