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

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

fast swap for BamRead?

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