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

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

refs #746. added function to modify one element in sequence

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 5.3 KB
Line 
1// $Id: BamRead.cc 2985 2013-02-18 07:36:14Z peter $
2//
3// Copyright (C) 2012, 2013 Peter Johansson
4//
5// This program is free software; you can redistribute it and/or modify
6// it under the terms of the GNU General Public License as published by
7// the Free Software Foundation; either version 3 of the License, or
8// (at your option) any later version.
9//
10// This program is distributed in the hope that it will be useful, but
11// WITHOUT ANY WARRANTY; without even the implied warranty of
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13// General Public License for more details.
14//
15// You should have received a copy of the GNU General Public License
16// along with this program. If not, see <http://www.gnu.org/licenses/>.
17
18#include <config.h>
19
20#include "BamRead.h"
21#include "config_bam.h"
22
23#include YAT_BAM_HEADER
24
25#include <algorithm>
26#include <cassert>
27#include <sstream>
28
29namespace theplu {
30namespace yat {
31namespace omic {
32
33  BamRead::BamRead(void)
34    : bam_(bam_init1())
35  {}
36
37
38  BamRead::BamRead(const BamRead& other)
39    : bam_(bam_dup1(other.bam_))
40  {
41  }
42
43
44  BamRead::~BamRead(void)
45  {
46    bam_destroy1(bam_);
47  }
48
49
50  const BamRead& BamRead::operator=(const BamRead& rhs)
51  {
52    if (bam_!=rhs.bam_)
53      bam_copy1(bam_, rhs.bam_);
54    return *this;
55  }
56
57
58  const uint8_t* BamRead::aux(void) const
59  {
60    return bam1_aux(bam_);
61  }
62
63
64  const uint32_t* BamRead::cigar(void) const
65  {
66    return bam1_cigar(bam_);
67  }
68
69
70  uint32_t BamRead::cigar(size_t i) const
71  {
72    assert(i<core().n_cigar);
73    return bam1_cigar(bam_)[i];
74  }
75
76
77  uint32_t BamRead::cigar_op(size_t i) const
78  {
79    assert(i<core().n_cigar);
80    return bam_cigar_op(cigar(i));
81  }
82
83
84  uint32_t BamRead::cigar_oplen(size_t i) const
85  {
86    assert(i<core().n_cigar);
87    return bam_cigar_oplen(cigar(i));
88  }
89
90
91  std::string BamRead::cigar_str(void) const
92  {
93    std::ostringstream os;
94    std::string str = BAM_CIGAR_STR;
95    for (uint32_t i=0; i<core().n_cigar; ++i) {
96      os << bam_cigar_oplen(bam1_cigar(bam_)[i])
97         << str[bam_cigar_op(bam1_cigar(bam_)[i])];
98    }
99    return os.str();
100  }
101
102
103  void BamRead::cigar(const std::vector<uint32_t>& c)
104  {
105    int m_data = bam_->m_data;
106    int data_len = bam_->data_len + 4*c.size() - 4*core().n_cigar;
107    // double capacity if needed
108    if (m_data < data_len) {
109      m_data = data_len;
110      kroundup32(m_data);
111    }
112    uint8_t* data = (uint8_t*)calloc(m_data, 1);
113    /*
114      data comes as qname, cigar, rest
115     */
116    // first copy qname
117    memcpy(data, bam_->data, core().l_qname);
118    // copy new cigar
119    if (c.size())
120      // each cigar element is 4 bytes
121      memcpy(data + core().l_qname, &c[0], 4*c.size());
122    // copy remaining data
123    memcpy(data + core().l_qname + 4*c.size(), bam1_seq(bam_),
124           data_len - 4*c.size() - core().l_qname);
125    bam_->m_data = m_data;
126    bam_->data_len = data_len;
127    core().n_cigar = c.size();
128    free(bam_->data);
129    bam_->data = data;
130  }
131
132
133  const bam1_core_t& BamRead::core(void) const
134  { return bam_->core; }
135
136
137  bam1_core_t& BamRead::core(void)
138  { return bam_->core; }
139
140
141  const char* BamRead::name(void) const
142  { return bam1_qname(bam_); }
143
144
145  uint8_t BamRead::qual(size_t i) const
146  { return *(bam1_qual(bam_)+i); }
147
148
149  std::string BamRead::sequence(void) const
150  {
151    // FIXME bam_nt16_rev_table not available in old lbam
152    std::string result(sequence_length(), ' ');
153    for (size_t i=0; i<result.size(); ++i)
154      result[i] = bam_nt16_rev_table[sequence(i)];
155    return result;
156  }
157
158
159  void BamRead::sequence(size_t i, uint8_t seq)
160  {
161    assert(bam_);
162    assert(seq < 16);
163    bam1_seq_seti(bam1_seq(bam_), i, seq);
164  }
165
166
167  uint8_t BamRead::sequence(size_t index) const
168  {
169    assert(bam_);
170    return bam1_seqi(bam1_seq(bam_),index);
171  }
172
173
174  uint32_t BamRead::sequence_length(void) const
175  {
176    assert(bam_);
177    return bam_->core.l_qseq;
178  }
179
180
181  int32_t BamRead::pos(void) const
182  {
183    assert(bam_);
184    return bam_->core.pos;
185  }
186
187
188  int32_t BamRead::tid(void) const
189  {
190    assert(bam_);
191    return bam_->core.tid;
192  }
193
194
195
196  int32_t BamRead::mpos(void) const
197  {
198    assert(bam_);
199    return bam_->core.mpos;
200  }
201
202
203  int32_t BamRead::end(void) const
204  {
205    assert(bam_);
206    return bam_calend(&core(), bam1_cigar(bam_));
207  }
208
209
210  void BamRead::swap(BamRead& other)
211  {
212    std::swap(bam_, other.bam_);
213  }
214
215
216  void swap(BamRead& lhs, BamRead& rhs)
217  {
218    lhs.swap(rhs);
219  }
220
221  bool same_query_name(const BamRead& lhs, const BamRead& rhs)
222  {
223    if (lhs.core().l_qname != rhs.core().l_qname)
224      return false;
225
226    return std::equal(lhs.name(),
227                      lhs.name() + lhs.core().l_qname,
228                      rhs.name());
229  }
230
231
232  bool soft_clipped(const BamRead& bam)
233  {
234    return left_soft_clipped(bam) || right_soft_clipped(bam);
235  }
236
237
238  uint32_t left_soft_clipped(const BamRead& bam)
239  {
240    for (size_t i=0; i<bam.core().n_cigar; ++i) {
241      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
242        return bam.cigar_oplen(i);
243      if (bam_cigar_type(bam.cigar_op(i))&1)
244        return 0;
245    }
246    return 0;
247  }
248
249
250  uint32_t right_soft_clipped(const BamRead& bam)
251  {
252    if (bam.core().n_cigar==0)
253      return 0;
254    // right clip can't be at first (most left cigar element)
255    for (size_t i=bam.core().n_cigar-1; i>0; --i) {
256      if (bam.cigar_op(i) == BAM_CSOFT_CLIP)
257        return bam.cigar_oplen(i);
258      if (bam_cigar_type(bam.cigar_op(i))&1)
259        return 0;
260    }
261    return 0;
262  }
263
264
265  bool BamLessPos::operator()(const BamRead& lhs, const BamRead& rhs) const
266  {
267    return lhs.tid() < rhs.tid() ||
268      (lhs.tid() == rhs.tid() && lhs.pos() < rhs.pos());
269  }
270
271
272  bool BamLessEnd::operator()(const BamRead& lhs, const BamRead& rhs) const
273  {
274    return lhs.tid() < rhs.tid() ||
275      (lhs.tid() == rhs.tid() && lhs.end() < rhs.end());
276  }
277
278}}}
Note: See TracBrowser for help on using the repository browser.