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

Last change on this file since 2972 was 2972, checked in by Peter, 10 years ago

refs #745. Implement a workaround when bam_nt16_rev_table is not
available in -lbam.

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