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

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

grammar in comment

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.1 KB
Line 
1// $Id: BamRead.cc 3324 2014-10-06 22:56:16Z 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 <cstring>
34#include <sstream>
35
36// bam1_seq_seti is not defined in 0.1.18 (or older), so
37// backport. Code is taken bam.h in version 0.1.19-dev.
38#ifndef bam1_seq_seti
39#define bam1_seq_seti(s, i, c) ( (s)[(i)>>1] = \
40  ((s)[(i)>>1] & 0xf<<(((i)&1)<<2)) | (c)<<((~(i)&1)<<2) )
41#endif
42
43namespace theplu {
44namespace yat {
45namespace omic {
46
47  BamRead::BamRead(void)
48    : bam_(bam_init1())
49  {}
50
51
52  BamRead::BamRead(const BamRead& other)
53    : bam_(bam_dup1(other.bam_))
54  {
55  }
56
57
58  BamRead::~BamRead(void)
59  {
60    bam_destroy1(bam_);
61  }
62
63
64  const BamRead& BamRead::operator=(const BamRead& rhs)
65  {
66    if (bam_!=rhs.bam_)
67      bam_copy1(bam_, rhs.bam_);
68    return *this;
69  }
70
71
72  const uint8_t* BamRead::aux(void) const
73  {
74    assert(bam_);
75    return bam1_aux(bam_);
76  }
77
78
79  const uint8_t* BamRead::aux(const char tag[2]) const
80  {
81    assert(bam_);
82    return bam_aux_get(bam_, tag);
83  }
84
85
86  void BamRead::aux_append(const char tag[2], char type,int len,uint8_t* data)
87  {
88    assert(bam_);
89    bam_aux_append(bam_, tag, type, len, data);
90  }
91
92
93  void BamRead::aux_del(const char tag[2])
94  {
95    assert(bam_);
96    uint8_t* s = bam_aux_get(bam_, tag);
97    if (!s) {
98      std::stringstream msg("BamRead::aux_del: absent tag: ");
99      msg << "'" << tag << "'";
100      throw utility::runtime_error(msg.str());
101    }
102    bam_aux_del(bam_, s);
103  }
104
105
106  int BamRead::aux_size(void) const
107  {
108    assert(bam_);
109    return bam_->l_aux;
110  }
111
112
113  uint32_t BamRead::calend(const bam1_core_t *c, const uint32_t *cigar) const
114  {
115#ifdef HAVE_BAM_CALEND
116    return bam_calend(&core(), bam1_cigar(bam_));
117#endif
118    int k, end = c->pos;
119    for (k = 0; k < c->n_cigar; ++k) {
120      int op  = bam_cigar_op(cigar[k]);
121      int len = bam_cigar_oplen(cigar[k]);
122      if (op == BAM_CBACK) { // move backward
123        int l, u, v;
124        if (k == c->n_cigar - 1) break; // skip trailing 'B'
125        for (l = k - 1, u = v = 0; l >= 0; --l) {
126          int op1  = bam_cigar_op(cigar[l]);
127          int len1 = bam_cigar_oplen(cigar[l]);
128          if (bam_cigar_type(op1)&1) { // consume query
129            if (u + len1 >= len) { // stop
130              if (bam_cigar_type(op1)&2) v += len - u;
131              break;
132            } else u += len1;
133          }
134          if (bam_cigar_type(op1)&2) v += len1;
135        }
136        end = l < 0? c->pos : end - v;
137      } else if (bam_cigar_type(op)&2) end += bam_cigar_oplen(cigar[k]);
138    }
139    return end;
140  }
141
142
143  const uint32_t* BamRead::cigar(void) const
144  {
145    return bam1_cigar(bam_);
146  }
147
148
149  uint32_t BamRead::cigar(size_t i) const
150  {
151    assert(i<core().n_cigar);
152    return bam1_cigar(bam_)[i];
153  }
154
155
156  uint32_t BamRead::cigar_op(size_t i) const
157  {
158    assert(i<core().n_cigar);
159    return bam_cigar_op(cigar(i));
160  }
161
162
163  uint32_t BamRead::cigar_oplen(size_t i) const
164  {
165    assert(i<core().n_cigar);
166    return bam_cigar_oplen(cigar(i));
167  }
168
169
170  std::string BamRead::cigar_str(void) const
171  {
172    std::ostringstream os;
173    for (uint32_t i=0; i<core().n_cigar; ++i) {
174      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
175         << bam_cigar_opchr(bam1_cigar(bam_)[i]);
176    }
177    return os.str();
178  }
179
180
181  void BamRead::cigar(const std::vector<uint32_t>& c)
182  {
183    // use the new function cigar(1) to implement this old deprecated variant
184    utility::Aligner::Cigar cig;
185    for (size_t i=0; i<c.size(); ++i)
186      cig.push_back(bam_cigar_op(c[i]), bam_cigar_oplen(c[i]));
187    cigar(cig);
188  }
189
190
191  void BamRead::cigar(const utility::Aligner::Cigar& c)
192  {
193    int offset = 4*c.size() - 4*core().n_cigar;
194    reserve(bam_->data_len + offset);
195    /*
196      data come as qname, cigar, rest
197     */
198    // first move remaining data to new location
199    if (offset)
200      memmove(bam1_seq(bam_)+offset, bam1_seq(bam_),
201              bam_->data_len - 4*core().n_cigar - core().l_qname);
202
203    // copy new cigar
204    for (size_t i=0; i<c.size(); ++i)
205      bam1_cigar(bam_)[i] = c[i];
206
207    bam_->data_len += offset;
208    core().n_cigar = c.size();
209    assert(bam_->data_len <= bam_->m_data);
210  }
211
212
213  const bam1_core_t& BamRead::core(void) const
214  { return bam_->core; }
215
216
217  bam1_core_t& BamRead::core(void)
218  { return bam_->core; }
219
220
221  int32_t BamRead::end(void) const
222  {
223    assert(bam_);
224    return this->calend(&core(), bam1_cigar(bam_));
225  }
226
227
228  uint16_t BamRead::flag(void) const
229  {
230    return bam_->core.flag;
231  }
232
233  const char* BamRead::name(void) const
234  { return bam1_qname(bam_); }
235
236
237  void BamRead::name(const std::string& n)
238  {
239    int offset = n.size() + 1 - core().l_qname;
240    assert(bam_);
241    reserve(bam_->data_len + offset);
242    // move remaining data
243    if (offset)
244      memmove(bam_->data + core().l_qname + offset,
245              bam_->data + core().l_qname,
246              bam_->data_len - core().l_qname);
247    core().l_qname += offset;
248    // copy new name
249    memcpy(bam1_qname(bam_), n.c_str(), n.size()+1);
250    bam_->data_len += offset;
251    assert(bam_->data_len <= bam_->m_data);
252  }
253
254
255  uint8_t BamRead::qual(size_t i) const
256  { return *(bam1_qual(bam_)+i); }
257
258
259  void BamRead::qual(size_t i, uint8_t q)
260  { *(bam1_qual(bam_)+i)=q; }
261
262
263  void BamRead::reserve(int n)
264  {
265    assert(bam_);
266    if (bam_->m_data >= n)
267      return;
268    // at least double capacity
269    bam_->m_data = std::max(bam_->m_data<<1, n);
270    bam_->data = (uint8_t*)realloc(bam_->data, bam_->m_data);
271    if (!bam_->data) {
272      throw utility::runtime_error("BamRead::reserve failed");
273    }
274  }
275
276
277  std::string BamRead::sequence(void) const
278  {
279    std::string result(sequence_length(), ' ');
280#ifdef HAVE_BAM_NT16_REV_TABLE
281    for (size_t i=0; i<result.size(); ++i)
282      result[i] = bam_nt16_rev_table[sequence(i)];
283#else
284    std::string table = "=ACMGRSVTWYHKDBN";
285    for (size_t i=0; i<result.size(); ++i)
286      result[i] = table[sequence(i)];
287#endif
288    return result;
289  }
290
291
292  void BamRead::sequence(const std::string& seq,
293                         const std::vector<uint8_t>& qual)
294  {
295    assert(seq.size()==qual.size());
296    // data consist of
297    // |-- qname -->-- cigar -->-- seq -->-- qual -->-- aux -->
298
299    int aux_pos = bam1_seq(bam_) - bam_->data + (seq.size()+1)/2 + qual.size();
300    int len = aux_pos + bam_->l_aux;
301    reserve(len);
302    // move aux to its new location
303    memmove(static_cast<void*>(bam_->data+aux_pos),
304            static_cast<const void*>(bam1_aux(bam_)), bam_->l_aux);
305    // copy sequence
306    for (size_t i=0; i<seq.size(); ++i)
307      sequence(i, bam_nt16_table[static_cast<size_t>(seq[i])]);
308    core().l_qseq = seq.size();
309    // copy quality
310    if (qual.size()) // FIXME is this needed or does memcpy work with n=0?
311      memcpy(static_cast<void*>(bam1_qual(bam_)),
312             static_cast<const void*>(&qual[0]), qual.size());
313    bam_->data_len = len;
314    assert(bam_->data_len <= bam_->m_data);
315  }
316
317
318  void BamRead::sequence(size_t i, uint8_t seq)
319  {
320    assert(bam_);
321    assert(seq < 16);
322    bam1_seq_seti(bam1_seq(bam_), i, seq);
323  }
324
325
326  uint8_t BamRead::sequence(size_t index) const
327  {
328    assert(bam_);
329    return bam1_seqi(bam1_seq(bam_),index);
330  }
331
332
333  uint32_t BamRead::sequence_length(void) const
334  {
335    assert(bam_);
336    return bam_->core.l_qseq;
337  }
338
339
340  int32_t BamRead::pos(void) const
341  {
342    assert(bam_);
343    return bam_->core.pos;
344  }
345
346
347  int32_t BamRead::tid(void) const
348  {
349    assert(bam_);
350    return bam_->core.tid;
351  }
352
353
354  int32_t BamRead::mtid(void) const
355  {
356    assert(bam_);
357    return bam_->core.mtid;
358  }
359
360
361  int32_t BamRead::mpos(void) const
362  {
363    assert(bam_);
364    return bam_->core.mpos;
365  }
366
367
368  void BamRead::swap(BamRead& other)
369  {
370    std::swap(bam_, other.bam_);
371  }
372
373
374  void swap(BamRead& lhs, BamRead& rhs)
375  {
376    lhs.swap(rhs);
377  }
378
379  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
380  {
381    if (lhs.core().l_qname != rhs.core().l_qname)
382      return false;
383
384    return !strcmp(lhs.name(), rhs.name());
385  }
386
387
388  bool soft_clipped(const BamRead& bam)
389  {
390    return left_soft_clipped(bam) || right_soft_clipped(bam);
391  }
392
393
394  uint32_t left_soft_clipped(const BamRead& bam)
395  {
396    for (size_t i=0; i<bam.core().n_cigar; ++i) {
397      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
398        return bam.cigar_oplen(i);
399      if (bam_cigar_type(bam.cigar_op(i))&1)
400        return 0;
401    }
402    return 0;
403  }
404
405
406  uint32_t right_soft_clipped(const BamRead& bam)
407  {
408    if (bam.core().n_cigar==0)
409      return 0;
410    // right clip can't be at first (most left cigar element)
411    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
412      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
413        return bam.cigar_oplen(i);
414      if (bam_cigar_type(bam.cigar_op(i))&1)
415        return 0;
416    }
417    return 0;
418  }
419
420
421  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
422  {
423    return lhs.tid() < rhs.tid() ||
424      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
425  }
426
427
428  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
429  {
430    return lhs.tid() < rhs.tid() ||
431      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
432  }
433
434}}}
Note: See TracBrowser for help on using the repository browser.