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

Last change on this file since 3167 was 3167, checked in by Peter, 8 years ago

avoid temporary variable

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 7.9 KB
Line 
1// $Id: BamRead.cc 3167 2014-01-21 07:56:03Z 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    for (uint32_t i=0; i<core().n_cigar; ++i) {
143      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
144         << bam_cigar_opchr(bam1_cigar(bam_)[i]);
145    }
146    return os.str();
147  }
148
149
150  void BamRead::cigar(const std::vector<uint32_t>& c)
151  {
152    int offset = 4*c.size() - 4*core().n_cigar;
153    reserve(bam_->data_len + offset);
154    /*
155      data comes as qname, cigar, rest
156     */
157    // first move remaining data to new location
158    if (offset)
159      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
160              bam_->data_len - 4*core().n_cigar - core().l_qname);
161
162    // copy new cigar
163    if (c.size())
164      // each cigar element is 4 bytes
165      memcpy(bam1_cigar(bam_), &c[0], 4*c.size());
166    bam_->data_len += offset;
167    core().n_cigar = c.size();
168    assert(bam_->data_len <= bam_->m_data);
169  }
170
171
172  const bam1_core_t& BamRead::core(void) const
173  { return bam_->core; }
174
175
176  bam1_core_t& BamRead::core(void)
177  { return bam_->core; }
178
179
180  const char* BamRead::name(void) const
181  { return bam1_qname(bam_); }
182
183
184  void BamRead::name(const std::string& n)
185  {
186    int offset = n.size() + 1 - core().l_qname;
187    assert(bam_);
188    reserve(bam_->data_len + offset);
189    // move remaining data
190    if (offset)
191      memmove(bam_->data + core().l_qname + offset,
192              bam_->data + core().l_qname,
193              bam_->data_len - core().l_qname);
194    core().l_qname += offset;
195    // copy new name
196    memcpy(bam1_qname(bam_), n.c_str(), n.size()+1);
197    bam_->data_len += offset;
198    assert(bam_->data_len <= bam_->m_data);
199  }
200
201
202  uint8_t BamRead::qual(size_t i) const
203  { return *(bam1_qual(bam_)+i); }
204
205
206  void BamRead::qual(size_t i, uint8_t q)
207  { *(bam1_qual(bam_)+i)=q; }
208
209
210  void BamRead::reserve(int n)
211  {
212    assert(bam_);
213    if (bam_->m_data >= n)
214      return;
215    // at least double capacity
216    bam_->m_data = std::max(bam_->m_data<<1, n);
217    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
218    if (!bam_->data) {
219      throw utility::runtime_error("BamRead::reserve failed");
220    }
221  }
222
223
224  std::string BamRead::sequence(void) const
225  {
226    std::string result(sequence_length(), ' ');
227#ifdef HAVE_BAM_NT16_REV_TABLE
228    for (size_t i=0; i<result.size(); ++i)
229      result[i] = bam_nt16_rev_table[sequence(i)];
230#else
231    std::string table = "=ACMGRSVTWYHKDBN";
232    for (size_t i=0; i<result.size(); ++i)
233      result[i] = table[sequence(i)];
234#endif
235    return result;
236  }
237
238
239  void BamRead::sequence(const std::string& seq,
240                         const std::vector<uint8_t>& qual)
241  {
242    assert(seq.size()==qual.size());
243    // data consist of
244    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
245
246    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
247    int len = aux_pos + bam_->l_aux;
248    reserve(len);
249    // move aux to its new location
250    memmove(static_cast<void*>(bam_->data+aux_pos),
251            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
252    // copy sequence
253    for (size_t i=0; i<seq.size(); ++i)
254      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
255    core().l_qseq = seq.size();
256    // copy quality
257    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
258      memcpy(static_cast<void*>(bam1_qual(bam_)),
259             static_cast<const void*>(&qual[0]), qual.size());
260    bam_->data_len = len;
261    assert(bam_->data_len <= bam_->m_data);
262  }
263
264
265  void BamRead::sequence(size_t i, uint8_t seq)
266  {
267    assert(bam_);
268    assert(seq < 16);
269    bam1_seq_seti(bam1_seq(bam_), i, seq);
270  }
271
272
273  uint8_t BamRead::sequence(size_t index) const
274  {
275    assert(bam_);
276    return bam1_seqi(bam1_seq(bam_),index);
277  }
278
279
280  uint32_t BamRead::sequence_length(void) const
281  {
282    assert(bam_);
283    return bam_->core.l_qseq;
284  }
285
286
287  int32_t BamRead::pos(void) const
288  {
289    assert(bam_);
290    return bam_->core.pos;
291  }
292
293
294  int32_t BamRead::tid(void) const
295  {
296    assert(bam_);
297    return bam_->core.tid;
298  }
299
300
301  int32_t BamRead::mtid(void) const
302  {
303    assert(bam_);
304    return bam_->core.mtid;
305  }
306
307
308  int32_t BamRead::mpos(void) const
309  {
310    assert(bam_);
311    return bam_->core.mpos;
312  }
313
314
315  int32_t BamRead::end(void) const
316  {
317    assert(bam_);
318    return bam_calend(&core(), bam1_cigar(bam_));
319  }
320
321
322  void BamRead::swap(BamRead& other)
323  {
324    std::swap(bam_, other.bam_);
325  }
326
327
328  void swap(BamRead& lhs, BamRead& rhs)
329  {
330    lhs.swap(rhs);
331  }
332
333  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
334  {
335    if (lhs.core().l_qname != rhs.core().l_qname)
336      return false;
337
338    return std::equal(lhs.name(),
339                      lhs.name() + lhs.core().l_qname,
340                      rhs.name());
341  }
342
343
344  bool soft_clipped(const BamRead& bam)
345  {
346    return left_soft_clipped(bam) || right_soft_clipped(bam);
347  }
348
349
350  uint32_t left_soft_clipped(const BamRead& bam)
351  {
352    for (size_t i=0; i<bam.core().n_cigar; ++i) {
353      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
354        return bam.cigar_oplen(i);
355      if (bam_cigar_type(bam.cigar_op(i))&1)
356        return 0;
357    }
358    return 0;
359  }
360
361
362  uint32_t right_soft_clipped(const BamRead& bam)
363  {
364    if (bam.core().n_cigar==0)
365      return 0;
366    // right clip can't be at first (most left cigar element)
367    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
368      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
369        return bam.cigar_oplen(i);
370      if (bam_cigar_type(bam.cigar_op(i))&1)
371        return 0;
372    }
373    return 0;
374  }
375
376
377  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
378  {
379    return lhs.tid() < rhs.tid() ||
380      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
381  }
382
383
384  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
385  {
386    return lhs.tid() < rhs.tid() ||
387      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
388  }
389
390}}}
Note: See TracBrowser for help on using the repository browser.