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

Last change on this file since 3306 was 3306, checked in by Peter, 7 years ago

merge ptach release 0.12.1 into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1// $Id: BamRead.cc 3306 2014-08-21 04:37:21Z peter $
2
3/*
4  Copyright (C) 2012, 2013, 2014 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  uint32_t BamRead::calend(const bam1_core_t *c, const uint32_t *cigar) const
113  {
114#ifdef HAVE_BAM_CALEND
115    return bam_calend(&core(), bam1_cigar(bam_));
116#endif
117    int k, end = c->pos;
118    for (k = 0; k < c->n_cigar; ++k) {
119      int op  = bam_cigar_op(cigar[k]);
120      int len = bam_cigar_oplen(cigar[k]);
121      if (op == BAM_CBACK) { // move backward
122        int l, u, v;
123        if (k == c->n_cigar - 1) break; // skip trailing 'B'
124        for (l = k - 1, u = v = 0; l >= 0; --l) {
125          int op1  = bam_cigar_op(cigar[l]);
126          int len1 = bam_cigar_oplen(cigar[l]);
127          if (bam_cigar_type(op1)&1) { // consume query
128            if (u + len1 >= len) { // stop
129              if (bam_cigar_type(op1)&2) v += len - u;
130              break;
131            } else u += len1;
132          }
133          if (bam_cigar_type(op1)&2) v += len1;
134        }
135        end = l < 0? c->pos : end - v;
136      } else if (bam_cigar_type(op)&2) end += bam_cigar_oplen(cigar[k]);
137    }
138    return end;
139  }
140
141
142  const uint32_t* BamRead::cigar(void) const
143  {
144    return bam1_cigar(bam_);
145  }
146
147
148  uint32_t BamRead::cigar(size_t i) const
149  {
150    assert(i<core().n_cigar);
151    return bam1_cigar(bam_)[i];
152  }
153
154
155  uint32_t BamRead::cigar_op(size_t i) const
156  {
157    assert(i<core().n_cigar);
158    return bam_cigar_op(cigar(i));
159  }
160
161
162  uint32_t BamRead::cigar_oplen(size_t i) const
163  {
164    assert(i<core().n_cigar);
165    return bam_cigar_oplen(cigar(i));
166  }
167
168
169  std::string BamRead::cigar_str(void) const
170  {
171    std::ostringstream os;
172    for (uint32_t i=0; i<core().n_cigar; ++i) {
173      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
174         << bam_cigar_opchr(bam1_cigar(bam_)[i]);
175    }
176    return os.str();
177  }
178
179
180  void BamRead::cigar(const std::vector<uint32_t>& c)
181  {
182    // use the new function cigar(1) to implement this old deprecated variant
183    utility::Aligner::Cigar cig;
184    for (size_t i=0; i<c.size(); ++i)
185      cig.push_back(bam_cigar_op(c[i]), bam_cigar_oplen(c[i]));
186    cigar(cig);
187  }
188
189
190  void BamRead::cigar(const utility::Aligner::Cigar& c)
191  {
192    int offset = 4*c.size() - 4*core().n_cigar;
193    reserve(bam_->data_len + offset);
194    /*
195      data comes as qname, cigar, rest
196     */
197    // first move remaining data to new location
198    if (offset)
199      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
200              bam_->data_len - 4*core().n_cigar - core().l_qname);
201
202    // copy new cigar
203    for (size_t i=0; i<c.size(); ++i)
204      bam1_cigar(bam_)[i] = c[i];
205
206    bam_->data_len += offset;
207    core().n_cigar = c.size();
208    assert(bam_->data_len <= bam_->m_data);
209  }
210
211
212  const bam1_core_t& BamRead::core(void) const
213  { return bam_->core; }
214
215
216  bam1_core_t& BamRead::core(void)
217  { return bam_->core; }
218
219
220  int32_t BamRead::end(void) const
221  {
222    assert(bam_);
223    return this->calend(&core(), bam1_cigar(bam_));
224  }
225
226
227  uint16_t BamRead::flag(void) const
228  {
229    return bam_->core.flag;
230  }
231
232  const char* BamRead::name(void) const
233  { return bam1_qname(bam_); }
234
235
236  void BamRead::name(const std::string& n)
237  {
238    int offset = n.size() + 1 - core().l_qname;
239    assert(bam_);
240    reserve(bam_->data_len + offset);
241    // move remaining data
242    if (offset)
243      memmove(bam_->data + core().l_qname + offset,
244              bam_->data + core().l_qname,
245              bam_->data_len - core().l_qname);
246    core().l_qname += offset;
247    // copy new name
248    memcpy(bam1_qname(bam_), n.c_str(), n.size()+1);
249    bam_->data_len += offset;
250    assert(bam_->data_len <= bam_->m_data);
251  }
252
253
254  uint8_t BamRead::qual(size_t i) const
255  { return *(bam1_qual(bam_)+i); }
256
257
258  void BamRead::qual(size_t i, uint8_t q)
259  { *(bam1_qual(bam_)+i)=q; }
260
261
262  void BamRead::reserve(int n)
263  {
264    assert(bam_);
265    if (bam_->m_data >= n)
266      return;
267    // at least double capacity
268    bam_->m_data = std::max(bam_->m_data<<1, n);
269    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
270    if (!bam_->data) {
271      throw utility::runtime_error("BamRead::reserve failed");
272    }
273  }
274
275
276  std::string BamRead::sequence(void) const
277  {
278    std::string result(sequence_length(), ' ');
279#ifdef HAVE_BAM_NT16_REV_TABLE
280    for (size_t i=0; i<result.size(); ++i)
281      result[i] = bam_nt16_rev_table[sequence(i)];
282#else
283    std::string table = "=ACMGRSVTWYHKDBN";
284    for (size_t i=0; i<result.size(); ++i)
285      result[i] = table[sequence(i)];
286#endif
287    return result;
288  }
289
290
291  void BamRead::sequence(const std::string& seq,
292                         const std::vector<uint8_t>& qual)
293  {
294    assert(seq.size()==qual.size());
295    // data consist of
296    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
297
298    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
299    int len = aux_pos + bam_->l_aux;
300    reserve(len);
301    // move aux to its new location
302    memmove(static_cast<void*>(bam_->data+aux_pos),
303            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
304    // copy sequence
305    for (size_t i=0; i<seq.size(); ++i)
306      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
307    core().l_qseq = seq.size();
308    // copy quality
309    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
310      memcpy(static_cast<void*>(bam1_qual(bam_)),
311             static_cast<const void*>(&qual[0]), qual.size());
312    bam_->data_len = len;
313    assert(bam_->data_len <= bam_->m_data);
314  }
315
316
317  void BamRead::sequence(size_t i, uint8_t seq)
318  {
319    assert(bam_);
320    assert(seq < 16);
321    bam1_seq_seti(bam1_seq(bam_), i, seq);
322  }
323
324
325  uint8_t BamRead::sequence(size_t index) const
326  {
327    assert(bam_);
328    return bam1_seqi(bam1_seq(bam_),index);
329  }
330
331
332  uint32_t BamRead::sequence_length(void) const
333  {
334    assert(bam_);
335    return bam_->core.l_qseq;
336  }
337
338
339  int32_t BamRead::pos(void) const
340  {
341    assert(bam_);
342    return bam_->core.pos;
343  }
344
345
346  int32_t BamRead::tid(void) const
347  {
348    assert(bam_);
349    return bam_->core.tid;
350  }
351
352
353  int32_t BamRead::mtid(void) const
354  {
355    assert(bam_);
356    return bam_->core.mtid;
357  }
358
359
360  int32_t BamRead::mpos(void) const
361  {
362    assert(bam_);
363    return bam_->core.mpos;
364  }
365
366
367  void BamRead::swap(BamRead& other)
368  {
369    std::swap(bam_, other.bam_);
370  }
371
372
373  void swap(BamRead& lhs, BamRead& rhs)
374  {
375    lhs.swap(rhs);
376  }
377
378  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
379  {
380    if (lhs.core().l_qname != rhs.core().l_qname)
381      return false;
382
383    return std::equal(lhs.name(),
384                      lhs.name() + lhs.core().l_qname,
385                      rhs.name());
386  }
387
388
389  bool soft_clipped(const BamRead& bam)
390  {
391    return left_soft_clipped(bam) || right_soft_clipped(bam);
392  }
393
394
395  uint32_t left_soft_clipped(const BamRead& bam)
396  {
397    for (size_t i=0; i<bam.core().n_cigar; ++i) {
398      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
399        return bam.cigar_oplen(i);
400      if (bam_cigar_type(bam.cigar_op(i))&1)
401        return 0;
402    }
403    return 0;
404  }
405
406
407  uint32_t right_soft_clipped(const BamRead& bam)
408  {
409    if (bam.core().n_cigar==0)
410      return 0;
411    // right clip can't be at first (most left cigar element)
412    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
413      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
414        return bam.cigar_oplen(i);
415      if (bam_cigar_type(bam.cigar_op(i))&1)
416        return 0;
417    }
418    return 0;
419  }
420
421
422  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
423  {
424    return lhs.tid() < rhs.tid() ||
425      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
426  }
427
428
429  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
430  {
431    return lhs.tid() < rhs.tid() ||
432      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
433  }
434
435}}}
Note: See TracBrowser for help on using the repository browser.