File: BamReader.cpp

package info (click to toggle)
libseqlib 1.1.2+dfsg-1~bpo9+1
  • links: PTS, VCS
  • area: main
  • in suites: stretch-backports
  • size: 1,460 kB
  • sloc: cpp: 7,176; sh: 805; makefile: 60
file content (400 lines) | stat: -rw-r--r-- 11,144 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
#include "SeqLib/BamReader.h"


//#define DEBUG_WALKER 1

namespace SeqLib {

// set the bam region
bool _Bam::SetRegion(const GenomicRegion& gp) {

  // mark it "open" again, may be new reads here
  mark_for_closure = false;
    
  //HTS set region 
  if ( (fp->format.format == 4 || fp->format.format == 6) && !idx)  // BAM (4) or CRAM (6)
    idx = SharedIndex(sam_index_load(fp.get(), m_in.c_str()), idx_delete());
  
  if (!idx) {
    if (m_in != "-")
      std::cerr << "Failed to load index for " << m_in << ". Rebuild samtools index" << std::endl;
    else
      std::cerr << "Random access with SetRegion not available for STDIN reading (no index file)" << std::endl;
    return false;
  }
  
  if (gp.chr >= m_hdr.NumSequences()) {
    std::cerr << "Failed to set region on " << gp << ". Chr ID is bigger than n_targets=" << m_hdr.NumSequences() << std::endl;
    return false;
  }
  
  // should work for BAM or CRAM
  hts_itr = SeqPointer<hts_itr_t>(sam_itr_queryi(idx.get(), gp.chr, gp.pos1, gp.pos2), hts_itr_delete());
  
  if (!hts_itr) {
    std::cerr << "Error: Failed to set region: " << gp << std::endl; 
    return false;
  }
  
  return true;
}

void BamReader::Reset() {
  for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b) 
     b->second.reset();
  m_region = GRC();
}

  bool BamReader::Reset(const std::string& f) {
    
    // cant reset what we don't have
    if (!m_bams.count(f))
      return false;
    m_bams[f].reset();
    return true;
}

  bool BamReader::Close() {
    
    bool success = true;
  for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b) 
      success = success && b->second.close();
    return success;
  }

  bool BamReader::Close(const std::string& f) {
    
    // cant close what we don't have
    if (!m_bams.count(f)) 
      return false;

    return m_bams[f].close();
  }

  /*  SharedHTSFile BamReader::GetHTSFile () const {
    if (!m_bams.size())
      throw std::runtime_error("No BAMs have been opened yet");
    return m_bams.begin()->second.fp;
  }

  SharedHTSFile BamReader::GetHTSFile(const std::string& f) const {
    _BamMap::const_iterator ff = m_bams.find(f);
    if (ff == m_bams.end())
      throw std::runtime_error("File " + f + " has not been opened yet");
    return ff->second.fp;
  }
  

  bool BamReader::SetPreloadedIndex(SharedIndex& i) {
    if (!m_bams.size())
      return false;
    m_bams.begin()->second.set_index(i);
    return true;
  }

  bool BamReader::SetPreloadedIndex(const std::string& f, SharedIndex& i) {
    if (!m_bams.count(f))
      return false;
    m_bams[f].set_index(i);
    return true;
  }

  */

  bool BamReader::SetRegion(const GenomicRegion& g) {
    m_region.clear();
    m_region.add(g);
    
    bool success = true;
    if (m_region.size()) {
      for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b) {
	b->second.m_region = &m_region;
	b->second.m_region_idx = 0; // set to the begining
	success = success && b->second.SetRegion(m_region[0]);
    }
    return success;
  }

  return false;
  
}

  bool BamReader::SetMultipleRegions(const GRC& grc) 
{
  if (grc.size() == 0) {
    std::cerr << "Warning: Trying to set an empty bam region"  << std::endl;
    return false;
  }
  
  m_region = grc;

  // go through and start all the BAMs at the first region
  bool success = true;
  if (m_region.size()) {
    for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b) {
      b->second.m_region = &m_region;
      b->second.m_region_idx = 0; // set to the begining
      success = success && b->second.SetRegion(m_region[0]);
    }
    return success;
  }
  
  return false;
}

  bool BamReader::Open(const std::string& bam) {

    // dont open same bam twice
    if (m_bams.count(bam))
      return false;
    
    _Bam new_bam(bam);
    new_bam.m_region = &m_region;
    bool success = new_bam.open_BAM_for_reading();
    m_bams.insert(std::pair<std::string, _Bam>(bam, new_bam));
    return success;
  }

  bool BamReader::Open(const std::vector<std::string>& bams) {
    
    bool pass = true;
    for (std::vector<std::string>::const_iterator i = bams.begin(); i != bams.end(); ++i)
      pass = pass && Open(*i);
    return pass;
  }
  
BamReader::BamReader() {}

  std::string BamReader::HeaderConcat() const {
    std::stringstream ss;
    for (_BamMap::const_iterator i = m_bams.begin(); i != m_bams.end(); ++i) 
      ss << i->second.m_hdr.AsString();
    return ss.str();

  }

  BamHeader BamReader::Header() const { 
    if (m_bams.size()) 
      return m_bams.begin()->second.m_hdr; 
    return BamHeader(); 
  }

  bool _Bam::open_BAM_for_reading() {

    // HTS open the reader
    fp = SharedHTSFile(hts_open(m_in.c_str(), "r"), htsFile_delete()); 

    // open cram reference
    if (!m_cram_reference.empty()) {
      char * m_cram_reference_cstr = strdup(m_cram_reference.c_str());
      int ret = cram_load_reference(fp->fp.cram, m_cram_reference_cstr);
      free(m_cram_reference_cstr);
      if (ret < 0) 
	throw std::invalid_argument("Could not read reference genome " + m_cram_reference + " for CRAM opt");
    }
   
    // check if opening failed
    if (!fp) 
      return false; 
    
    // read the header and create a BamHeader
    bam_hdr_t * hdr = sam_hdr_read(fp.get());
    m_hdr = BamHeader(hdr); // calls BamHeader(bam_hdr_t), makes a copy

    // deallocate the memory we just made
    if (hdr)
      bam_hdr_destroy(hdr);
    
    // if BAM header opening failed, return false
    if (!m_hdr.get()) 
      return false;
    
    // everything worked
    return true;
    
  }

  void BamReader::SetCramReference(const std::string& ref) {
    m_cram_reference = ref;
    for (_BamMap::iterator b = m_bams.begin(); b != m_bams.end(); ++b)
      b->second.m_cram_reference = ref;
  }

bool BamReader::GetNextRecord(BamRecord& r) {

  // shortcut if we have only a single bam
  if (m_bams.size() == 1) {
    
    if (m_bams.begin()->second.fp.get() == NULL || m_bams.begin()->second.mark_for_closure) // cant read if not opened
      return false;
    
    // try and get the next read
    int32_t status = m_bams.begin()->second.load_read(r);
    if (status >= 0)
      return true;
    if (status == -1) {
      // didn't find anything, clear it
      m_bams.begin()->second.mark_for_closure = true;
      return false;
    }
    
    // run time error
    std::stringstream ss;
    ss << "sam_read1 return status: " << status << " file: " << m_bams.begin()->first;
    throw std::runtime_error(ss.str());
    return false;
  }

  // loop the files and load the next read
  // for the one that was emptied last
  for (_BamMap::iterator bam = m_bams.begin(); bam != m_bams.end(); ++bam) {

    _Bam *tb = &(bam->second);

    // if marked, then don't even try on this BAM
    if (tb->mark_for_closure)
      continue;

    // skip un-opened BAMs
    if (tb->fp.get() == NULL) 
      continue;

    // if next read is not marked as empty, skip to next
    if (!tb->empty)
      continue; 
    
    // load the next read
    int32_t status = tb->load_read(r);
    if (status == -1) {
      // can't load, so mark for closing
      tb->empty = true;
      tb->mark_for_closure = true; // no more reads in this BAM
      continue; 
    } else if (status < 0) { // error sent back from sam_read1
      // run time error
      std::stringstream ss;
      ss << "sam_read1 return status: " << status << " file: " << bam->first;
      throw std::runtime_error(ss.str());
    }
    
  }

  // for multiple bams, choose the one to return
  // sort based on chr and left-most alignment pos. Same as samtools
  int min_chr = INT_MAX;
  int min_pos = INT_MAX;
  _BamMap::iterator hit; 
  bool found = false; // did we find a valid read

  for (_BamMap::iterator bam = m_bams.begin(); 
       bam != m_bams.end(); ++bam) {

    // dont check if already marked for removal or doesnt need new read
    if (bam->second.empty || bam->second.mark_for_closure) 
      continue;

    found = true;
    if (bam->second.next_read.ChrID() < min_chr || // check if read in this BAM is lowest
	(bam->second.next_read.Position() <  min_pos && bam->second.next_read.ChrID() == min_chr)) {
      min_pos = bam->second.next_read.Position();
      min_chr = bam->second.next_read.ChrID();
      hit = bam; // read is lowest, so mark this BAM as having the hit
    }
  }

  // mark the one we just found as empty
  if (found) {
    r = hit->second.next_read; // read is lowest, so assign
    hit->second.empty = true;  // mark as empty, so we fill this slot again
  }
  
  return found;
}
  
std::string BamReader::PrintRegions() const {

  std::stringstream ss;
  //for (GRC::const_iterator r = m_region.begin(); r != m_region.end(); ++r)
  //  ss << *r << std::endl;
  return(ss.str());

}

  int32_t _Bam::load_read(BamRecord& r) {

  // allocated the memory
  bam1_t* b = bam_init1(); 
  int32_t valid = -1; // start with EOF return code

  if (hts_itr.get() == NULL) {
    valid = sam_read1(fp.get(), m_hdr.get_(), b);    

    if (valid < 0) { 
      
#ifdef DEBUG_WALKER
      std::cerr << "ended reading on null hts_itr" << std::endl;
#endif
      //goto endloop;
      bam_destroy1(b);
      return valid;
    }
  } else {
    
    //changed to sam from hts_itr_next
    // move to next region of bam
    valid = sam_itr_next(fp.get(), hts_itr.get(), b);
  }
  
  if (valid < 0) { // read still not found
    do {
      
#ifdef DEBUG_WALKER
      std::cerr << "Failed read, trying next region. Moving counter to " << m_region_idx << " of " << m_region.size() << " FP: "  << fp_htsfile << " hts_itr " << std::endl;
#endif
      // try next region, return if no others to try
      ++m_region_idx; // increment to next region
      if (m_region_idx >= m_region->size()) {
	bam_destroy1(b);
	return valid;
      }
	//goto endloop;
      
      // next region exists, try it
      SetRegion(m_region->at(m_region_idx));
      valid = sam_itr_next(fp.get(), hts_itr.get(), b);
    } while (valid <= 0); // keep trying regions until works
  }
  
  // if we got here, then we found a read in this BAM
  empty = false;
  next_read.assign(b); // assign the shared_ptr for the bam1_t
  r = next_read;

  return valid;
}

std::ostream& operator<<(std::ostream& out, const BamReader& b)
{
  for(_BamMap::const_iterator bam = b.m_bams.begin(); bam != b.m_bams.end(); ++bam)
    out << ":" << bam->second.GetFileName() << std::endl; 

  if (b.m_region.size() && b.m_region.size() < 20) {
    out << " ------- BamReader Regions ----------" << std::endl;;
    //for (GRC::const_iterator r = b.m_region.begin(); r != b.m_region.end(); ++r)
    //  out << *i << std::endl;
  } 
  else if (b.m_region.size() >= 20) {
    int wid = 0;
    //for (GRC::const_iterator r = b.m_region.begin(); r != b.m_region.end(); ++r)
    //  wid += r->Width();
    out << " ------- BamReader Regions ----------" << std::endl;;
    out << " -- " << b.m_region.size() << " regions covering " << AddCommas(wid) << " bp of sequence"  << std::endl;
  }
  else 
    out << " - BamReader - Walking whole genome -" << std::endl;

  out <<   " ------------------------------------";
  return out;
}

}