source: trunk/yat/omic/BamRead.cc @ 2986

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

merge 0.10-stable branch into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 5.5 KB
Line 
1// $Id: BamRead.cc 2986 2013-02-18 08:07:44Z peter $
2//
3// Copyright (C) 2012, 2013 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    std::string result(sequence_length(), ' ');
152#ifdef HAVE_BAM_NT16_REV_TABLE
153    for (size_t i=0; i<result.size(); ++i)
154      result[i] = bam_nt16_rev_table[sequence(i)];
155#else
156    std::string table = "=ACMGRSVTWYHKDBN";
157    for (size_t i=0; i<result.size(); ++i)
158      result[i] = table[sequence(i)];
159#endif
160    return result;
161  }
162
163
164  void BamRead::sequence(size_t i, uint8_t seq)
165  {
166    assert(bam_);
167    assert(seq < 16);
168    bam1_seq_seti(bam1_seq(bam_), i, seq);
169  }
170
171
172  uint8_t BamRead::sequence(size_t index) const
173  {
174    assert(bam_);
175    return bam1_seqi(bam1_seq(bam_),index);
176  }
177
178
179  uint32_t BamRead::sequence_length(void) const
180  {
181    assert(bam_);
182    return bam_->core.l_qseq;
183  }
184
185
186  int32_t BamRead::pos(void) const
187  {
188    assert(bam_);
189    return bam_->core.pos;
190  }
191
192
193  int32_t BamRead::tid(void) const
194  {
195    assert(bam_);
196    return bam_->core.tid;
197  }
198
199
200  int32_t BamRead::mtid(void) const
201  {
202    assert(bam_);
203    return bam_->core.mtid;
204  }
205
206
207  int32_t BamRead::mpos(void) const
208  {
209    assert(bam_);
210    return bam_->core.mpos;
211  }
212
213
214  int32_t BamRead::end(void) const
215  {
216    assert(bam_);
217    return bam_calend(&core(), bam1_cigar(bam_));
218  }
219
220
221  void BamRead::swap(BamRead& other)
222  {
223    std::swap(bam_, other.bam_);
224  }
225
226
227  void swap(BamRead& lhs, BamRead& rhs)
228  {
229    lhs.swap(rhs);
230  }
231
232  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
233  {
234    if (lhs.core().l_qname != rhs.core().l_qname)
235      return false;
236
237    return std::equal(lhs.name(),
238                      lhs.name() + lhs.core().l_qname,
239                      rhs.name());
240  }
241
242
243  bool soft_clipped(const BamRead& bam)
244  {
245    return left_soft_clipped(bam) || right_soft_clipped(bam);
246  }
247
248
249  uint32_t left_soft_clipped(const BamRead& bam)
250  {
251    for (size_t i=0; i<bam.core().n_cigar; ++i) {
252      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
253        return bam.cigar_oplen(i);
254      if (bam_cigar_type(bam.cigar_op(i))&1)
255        return 0;
256    }
257    return 0;
258  }
259
260
261  uint32_t right_soft_clipped(const BamRead& bam)
262  {
263    if (bam.core().n_cigar==0)
264      return 0;
265    // right clip can't be at first (most left cigar element)
266    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
267      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
268        return bam.cigar_oplen(i);
269      if (bam_cigar_type(bam.cigar_op(i))&1)
270        return 0;
271    }
272    return 0;
273  }
274
275
276  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
277  {
278    return lhs.tid() < rhs.tid() ||
279      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
280  }
281
282
283  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
284  {
285    return lhs.tid() < rhs.tid() ||
286      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
287  }
288
289}}}
Note: See TracBrowser for help on using the repository browser.