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

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

closes #746. New function BamRead::name(1)

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