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

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

new function BamRead::sequence(2). refs #746

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.8 KB
Line 
1// $Id: BamRead.cc 3031 2013-04-27 11:24:13Z peter $
2
3/*
4  Copyright (C) 2012, 2013 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with this program. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <config.h>
23
24#include "BamRead.h"
25#include "config_bam.h"
26
27#include YAT_BAM_HEADER
28
29#include "yat/utility/Exception.h"
30
31#include <algorithm>
32#include <cassert>
33#include <sstream>
34
35// bam1_seq_seti is not defined in 0.1.18 (or older), so
36// backport. Code is taken bam.h in version 0.1.19-dev.
37#ifndef bam1_seq_seti
38#define bam1_seq_seti(s, i, c) ( (s)[(i)>>1] = \
39  ((s)[(i)>>1] & 0xf<<(((i)&1)<<2)) | (c)<<((~(i)&1)<<2) )
40#endif
41
42namespace theplu {
43namespace yat {
44namespace omic {
45
46  BamRead::BamRead(void)
47    : bam_(bam_init1())
48  {}
49
50
51  BamRead::BamRead(const BamRead& other)
52    : bam_(bam_dup1(other.bam_))
53  {
54  }
55
56
57  BamRead::~BamRead(void)
58  {
59    bam_destroy1(bam_);
60  }
61
62
63  const BamRead& BamRead::operator=(const BamRead& rhs)
64  {
65    if (bam_!=rhs.bam_)
66      bam_copy1(bam_, rhs.bam_);
67    return *this;
68  }
69
70
71  const uint8_t* BamRead::aux(void) const
72  {
73    assert(bam_);
74    return bam1_aux(bam_);
75  }
76
77
78  const uint8_t* BamRead::aux(const char tag[2]) const
79  {
80    assert(bam_);
81    return bam_aux_get(bam_, tag);
82  }
83
84
85  void BamRead::aux_append(const char tag[2], char type,int len,uint8_t* data)
86  {
87    assert(bam_);
88    bam_aux_append(bam_, tag, type, len, data);
89  }
90
91
92  void BamRead::aux_del(const char tag[2])
93  {
94    assert(bam_);
95    uint8_t* s = bam_aux_get(bam_, tag);
96    if (!s) {
97      std::stringstream msg("BamRead::aux_del: absent tag: ");
98      msg << "'" << tag << "'";
99      throw utility::runtime_error(msg.str());
100    }
101    bam_aux_del(bam_, s);
102  }
103
104
105  int BamRead::aux_size(void) const
106  {
107    assert(bam_);
108    return bam_->l_aux;
109  }
110
111
112  const uint32_t* BamRead::cigar(void) const
113  {
114    return bam1_cigar(bam_);
115  }
116
117
118  uint32_t BamRead::cigar(size_t i) const
119  {
120    assert(i<core().n_cigar);
121    return bam1_cigar(bam_)[i];
122  }
123
124
125  uint32_t BamRead::cigar_op(size_t i) const
126  {
127    assert(i<core().n_cigar);
128    return bam_cigar_op(cigar(i));
129  }
130
131
132  uint32_t BamRead::cigar_oplen(size_t i) const
133  {
134    assert(i<core().n_cigar);
135    return bam_cigar_oplen(cigar(i));
136  }
137
138
139  std::string BamRead::cigar_str(void) const
140  {
141    std::ostringstream os;
142    std::string str = BAM_CIGAR_STR;
143    for (uint32_t i=0; i<core().n_cigar; ++i) {
144      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
145         << str[bam_cigar_op(bam1_cigar(bam_)[i])];
146    }
147    return os.str();
148  }
149
150
151  void BamRead::cigar(const std::vector<uint32_t>& c)
152  {
153    // special case when cigar size is not changed
154    if (c.size() == core().n_cigar) {
155      if (c.size())
156        memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
157      return;
158    }
159
160    int m_data = bam_->m_data;
161    int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
162    // double capacity if needed
163    if (m_data < data_len) {
164      m_data = data_len;
165      kroundup32(m_data);
166    }
167    uint8_t* data = (uint8_t*)calloc(m_data, 1);
168    /*
169      data comes as qname, cigar, rest
170     */
171    // first copy qname
172    memcpy(data, bam_->data, core().l_qname);
173    // copy new cigar
174    if (c.size())
175      // each cigar element is 4 bytes
176      memcpy(data + core().l_qname, &c[0], 4*c.size());
177    // copy remaining data
178    memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
179           data_len - 4*c.size() - core().l_qname);
180    bam_->m_data = m_data;
181    bam_->data_len = data_len;
182    core().n_cigar = c.size();
183    free(bam_->data);
184    bam_->data = data;
185  }
186
187
188  const bam1_core_t& BamRead::core(void) const
189  { return bam_->core; }
190
191
192  bam1_core_t& BamRead::core(void)
193  { return bam_->core; }
194
195
196  const char* BamRead::name(void) const
197  { return bam1_qname(bam_); }
198
199
200  uint8_t BamRead::qual(size_t i) const
201  { return *(bam1_qual(bam_)+i); }
202
203
204  void BamRead::qual(size_t i, uint8_t q)
205  { *(bam1_qual(bam_)+i)=q; }
206
207
208  void BamRead::reserve(int n)
209  {
210    assert(bam_);
211    if (bam_->m_data >= n)
212      return;
213    // at least double capacity
214    bam_->m_data = std::max(bam_->m_data<<1, n);
215    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
216    if (!bam_->data) {
217      throw utility::runtime_error("BamRead::reserve failed");
218    }
219  }
220
221
222  std::string BamRead::sequence(void) const
223  {
224    std::string result(sequence_length(), ' ');
225#ifdef HAVE_BAM_NT16_REV_TABLE
226    for (size_t i=0; i<result.size(); ++i)
227      result[i] = bam_nt16_rev_table[sequence(i)];
228#else
229    std::string table = "=ACMGRSVTWYHKDBN";
230    for (size_t i=0; i<result.size(); ++i)
231      result[i] = table[sequence(i)];
232#endif
233    return result;
234  }
235
236
237  void BamRead::sequence(const std::string& seq,
238                         const std::vector<uint8_t>& qual)
239  {
240    assert(seq.size()==qual.size());
241    // data consist of
242    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
243
244    uint8_t* new_aux = bam1_seq(bam_) + (seq.size()+1)/2 + qual.size();
245
246    int len = new_aux + bam_->l_aux - bam_->data;
247    reserve(len);
248    // move aux to its new location
249    memmove(static_cast<void*>(new_aux),
250            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
251    // copy sequence
252    for (size_t i=0; i<seq.size(); ++i)
253      sequence(i, bam_nt16_table[seq[i]]);
254    core().l_qseq = seq.size();
255    // copy quality
256    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
257      memcpy(static_cast<void*>(bam1_qual(bam_)),
258             static_cast<const void*>(&qual[0]), qual.size());
259  }
260
261
262  void BamRead::sequence(size_t i, uint8_t seq)
263  {
264    assert(bam_);
265    assert(seq < 16);
266    bam1_seq_seti(bam1_seq(bam_), i, seq);
267  }
268
269
270  uint8_t BamRead::sequence(size_t index) const
271  {
272    assert(bam_);
273    return bam1_seqi(bam1_seq(bam_),index);
274  }
275
276
277  uint32_t BamRead::sequence_length(void) const
278  {
279    assert(bam_);
280    return bam_->core.l_qseq;
281  }
282
283
284  int32_t BamRead::pos(void) const
285  {
286    assert(bam_);
287    return bam_->core.pos;
288  }
289
290
291  int32_t BamRead::tid(void) const
292  {
293    assert(bam_);
294    return bam_->core.tid;
295  }
296
297
298  int32_t BamRead::mtid(void) const
299  {
300    assert(bam_);
301    return bam_->core.mtid;
302  }
303
304
305  int32_t BamRead::mpos(void) const
306  {
307    assert(bam_);
308    return bam_->core.mpos;
309  }
310
311
312  int32_t BamRead::end(void) const
313  {
314    assert(bam_);
315    return bam_calend(&core(), bam1_cigar(bam_));
316  }
317
318
319  void BamRead::swap(BamRead& other)
320  {
321    std::swap(bam_, other.bam_);
322  }
323
324
325  void swap(BamRead& lhs, BamRead& rhs)
326  {
327    lhs.swap(rhs);
328  }
329
330  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
331  {
332    if (lhs.core().l_qname != rhs.core().l_qname)
333      return false;
334
335    return std::equal(lhs.name(),
336                      lhs.name() + lhs.core().l_qname,
337                      rhs.name());
338  }
339
340
341  bool soft_clipped(const BamRead& bam)
342  {
343    return left_soft_clipped(bam) || right_soft_clipped(bam);
344  }
345
346
347  uint32_t left_soft_clipped(const BamRead& bam)
348  {
349    for (size_t i=0; i<bam.core().n_cigar; ++i) {
350      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
351        return bam.cigar_oplen(i);
352      if (bam_cigar_type(bam.cigar_op(i))&1)
353        return 0;
354    }
355    return 0;
356  }
357
358
359  uint32_t right_soft_clipped(const BamRead& bam)
360  {
361    if (bam.core().n_cigar==0)
362      return 0;
363    // right clip can't be at first (most left cigar element)
364    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
365      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
366        return bam.cigar_oplen(i);
367      if (bam_cigar_type(bam.cigar_op(i))&1)
368        return 0;
369    }
370    return 0;
371  }
372
373
374  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
375  {
376    return lhs.tid() < rhs.tid() ||
377      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
378  }
379
380
381  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
382  {
383    return lhs.tid() < rhs.tid() ||
384      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
385  }
386
387}}}
Note: See TracBrowser for help on using the repository browser.