Class: HTS::Bam

Inherits:
Hts
  • Object
show all
Includes:
Enumerable
Defined in:
lib/hts/bam.rb,
lib/hts/bam/auxi.rb,
lib/hts/bam/flag.rb,
lib/hts/bam/cigar.rb,
lib/hts/bam/header.rb,
lib/hts/bam/pileup.rb,
lib/hts/bam/record.rb,
lib/hts/bam/mpileup.rb,
lib/hts/bam/base_mod.rb,
lib/hts/bam/header_record.rb

Overview

A class for working with SAM, BAM, CRAM files.

Defined Under Namespace

Classes: Aux, BaseMod, Cigar, Flag, Header, HeaderRecord, MissingIndexError, Mpileup, OpenError, Pileup, ReadError, Record, WriteError

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(file_name, mode = "r", index: nil, fai: nil, threads: nil, build_index: false) ⇒ Bam

Returns a new instance of Bam.



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
 
# File 'lib/hts/bam.rb', line 72

def initialize(file_name, mode = "r", index: nil, fai: nil, threads: nil,
               build_index: false)
  if block_given?
    message = "HTS::Bam.new() does not take block; Please use HTS::Bam.open() instead"
    raise message
  end

  # NOTE: Do not check for the existence of local files, since file_names may be remote URIs.

  @file_name  = file_name
  @index_name = index
  @mode       = mode
  @nthreads   = threads
  @index_load_attempted = false
  @native = Native::BamFileHandle.open(@file_name, mode)

  # Auto-detect and set reference for CRAM files
  if fai.nil? && @file_name.end_with?(".cram")
    # Try to find reference file in the same directory
    base_name = File.basename(@file_name, ".cram")
    dir_name = File.dirname(@file_name)
    potential_ref = File.join(dir_name, "#{base_name}.fa")

    # For remote URLs, assume reference exists; for local files, check existence
    fai = potential_ref if @file_name.start_with?("http") || File.exist?(potential_ref)
  end

  if fai
    r = @native.set_fai(fai)
    raise "Failed to load fasta index: #{fai}" if r < 0
  end

  set_threads(threads) if threads

  if writing?
    @auto_index_on_close = build_index
    @index_name_on_close = index
    return
  end

  @header = Bam::Header.new(@native.read_header)
  if build_index
    build_index(index)
    load_index(index)
  elsif index
    load_index(index)
  end
  @start_position = tell
end

Instance Attribute Details

#file_name ⇒ Object (readonly)

Returns the value of attribute file_name.



39
40
41
 
# File 'lib/hts/bam.rb', line 39

def file_name
  @file_name
end

#header ⇒ Object

Returns the value of attribute header.



39
40
41
 
# File 'lib/hts/bam.rb', line 39

def header
  @header
end

#index_name ⇒ Object (readonly)

Returns the value of attribute index_name.



39
40
41
 
# File 'lib/hts/bam.rb', line 39

def index_name
  @index_name
end

#mode ⇒ Object (readonly)

Returns the value of attribute mode.



39
40
41
 
# File 'lib/hts/bam.rb', line 39

def mode
  @mode
end

#nthreads ⇒ Object (readonly)

Returns the value of attribute nthreads.



39
40
41
 
# File 'lib/hts/bam.rb', line 39

def nthreads
  @nthreads
end

Class Method Details

.build_index(file_name, index_name = nil, min_shift = 0, threads = 0, verbose = true) ⇒ Object 



53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
 
# File 'lib/hts/bam.rb', line 53

def self.build_index(file_name, index_name = nil, min_shift = 0, threads = 0, verbose = true)
  if verbose
    if index_name
      warn "Create index for #{file_name} to #{index_name}"
    else
      warn "Create index for #{file_name}"
    end
  end

  case Native::BamFileHandle.build_index(file_name, index_name, min_shift, threads)
  when 0 # successful
  when -1 then raise "indexing failed"
  when -2 then raise "opening #{file_name} failed"
  when -3 then raise "format not indexable"
  when -4 then raise "failed to create and/or save the index"
  else raise "unknown error"
  end
end

.filter_records(records, required_flags: 0, excluded_flags: 0, min_mapq: 0, tid: nil, beg: nil, end_: nil) ⇒ Object

Filter an owning batch of records in one native pass when available.



26
27
28
29
30
31
32
33
34
35
36
37
 
# File 'lib/hts/bam.rb', line 26

def self.filter_records(records, required_flags: 0, excluded_flags: 0,
                        min_mapq: 0, tid: nil, beg: nil, end_: nil)
  required_flags = Integer(required_flags)
  excluded_flags = Integer(excluded_flags)
  min_mapq = Integer(min_mapq)
  Array(records).select do |record|
    flags = record.flag_value
    (flags & required_flags) == required_flags && (flags & excluded_flags).zero? &&
      record.mapq >= min_mapq && (tid.nil? || record.tid == Integer(tid)) &&
      (beg.nil? || record.endpos > Integer(beg)) && (end_.nil? || record.pos < Integer(end_))
  end
end

.open(*args, **kw) ⇒ Object 



41
42
43
44
45
46
47
48
49
50
51
 
# File 'lib/hts/bam.rb', line 41

def self.open(*args, **kw)
  file = new(*args, **kw) # do not yield
  return file unless block_given?

  begin
    result = yield file
  ensure
    file.close
  end
  result
end

Instance Method Details

#<<(record) ⇒ Object 



222
223
224
225
 
# File 'lib/hts/bam.rb', line 222

def <<(record)
  write(record)
  self
end

#aux(tag) ⇒ Object Also known as: aux_array

FIXME: experimental



247
248
249
250
251
252
253
254
 
# File 'lib/hts/bam.rb', line 247

def aux(tag)
  check_closed

  position = tell
  ary = map { |r| r.aux(tag) }
  seek(position) if position
  ary
end

#build_index(index_name = nil, min_shift: 0, verbose: true) ⇒ Object 



122
123
124
125
126
127
128
129
 
# File 'lib/hts/bam.rb', line 122

def build_index(index_name = nil, min_shift: 0, verbose: true)
  check_closed

  self.class.build_index(@file_name, index_name, min_shift, @nthreads || 0, verbose)
  @index_name = index_name
  @index_load_attempted = false
  self # for method chaining
end

#chrom ⇒ Array

Get chrom array

Returns:

  • (Array) —

    the chrom array



233
 
# File 'lib/hts/bam.rb', line 233

define_getter :chrom

#cigar ⇒ Array

Get cigar array

Returns:

  • (Array) —

    the cigar array



236
 
# File 'lib/hts/bam.rb', line 236

define_getter :cigar

#close ⇒ Object

Raises:



151
152
153
154
155
156
157
158
159
160
161
162
 
# File 'lib/hts/bam.rb', line 151

def close
  was_closed = closed?
  result = @native&.close
  raise WriteError, "Failed to close #{@file_name}: buffered output may be incomplete" if writing? && result&.negative?

  if writing? && @auto_index_on_close && !was_closed
    @auto_index_on_close = false
    self.class.build_index(@file_name, @index_name_on_close, 0, @nthreads || 0, false)
  end

  nil
end

#closed? ⇒ Boolean

Returns:

  • (Boolean) 


164
 
# File 'lib/hts/bam.rb', line 164

def closed? = @native.nil? || @native.closed?

#collect_records ⇒ Object

Materialize independent records from the current stream position. Unlike each.to_a, every element owns its bam1_t storage.



259
260
261
 
# File 'lib/hts/bam.rb', line 259

def collect_records
  each(copy: true).to_a
end

#each(copy: false, &block) ⇒ Object

Iterate alignment records in this file.

Performance and memory semantics:

  • copy: false (default) reuses a single Record instance and its underlying bam1_t buffer. The yielded Record MUST NOT be stored beyond the block; its content will be overwritten by the next iteration. If you need to retain it, call rec = rec.dup.
  • copy: true yields a fresh Record per iteration (deep-copied via bam_dup1). Slower, safe to keep.


307
308
309
310
311
312
313
 
# File 'lib/hts/bam.rb', line 307

def each(copy: false, &block)
  if copy
    each_record_copy(&block)
  else
    each_record_reuse(&block)
  end
end

#each_aux(tag) ⇒ Object

FIXME: experimental



289
290
291
292
293
294
295
296
297
298
 
# File 'lib/hts/bam.rb', line 289

def each_aux(tag)
  check_closed
  return to_enum(__method__, tag) unless block_given?

  each do |record|
    yield record.aux(tag)
  end

  self
end

#each_chrom ⇒ Object

Get chrom iterator



275
 
# File 'lib/hts/bam.rb', line 275

define_iterator :chrom

#each_cigar ⇒ Object

Get cigar iterator



278
 
# File 'lib/hts/bam.rb', line 278

define_iterator :cigar

#each_flag ⇒ Object

Get flag iterator



274
 
# File 'lib/hts/bam.rb', line 274

define_iterator :flag

#each_insert_size ⇒ Object Also known as: each_isize

Get insert_size iterator



281
 
# File 'lib/hts/bam.rb', line 281

define_iterator :insert_size

#each_mapq ⇒ Object

Get mapq iterator



277
 
# File 'lib/hts/bam.rb', line 277

define_iterator :mapq

#each_mate_chrom ⇒ Object

Get mate_chrom iterator



279
 
# File 'lib/hts/bam.rb', line 279

define_iterator :mate_chrom

#each_mate_pos ⇒ Object Also known as: each_mpos

Get mate_pos iterator



280
 
# File 'lib/hts/bam.rb', line 280

define_iterator :mate_pos

#each_pos ⇒ Object

Get pos iterator



276
 
# File 'lib/hts/bam.rb', line 276

define_iterator :pos

#each_qname ⇒ Object

Get qname iterator



273
 
# File 'lib/hts/bam.rb', line 273

define_iterator :qname

#each_qual ⇒ Object

Get qual iterator



283
 
# File 'lib/hts/bam.rb', line 283

define_iterator :qual

#each_seq ⇒ Object

Get seq iterator



282
 
# File 'lib/hts/bam.rb', line 282

define_iterator :seq

#file_format ⇒ Object 



165
 
# File 'lib/hts/bam.rb', line 165

def file_format = @native.file_format

#file_format_version ⇒ Object 



166
 
# File 'lib/hts/bam.rb', line 166

def file_format_version = @native.file_format_version

#flag ⇒ Array

Get flag array

Returns:

  • (Array) —

    the flag array



232
 
# File 'lib/hts/bam.rb', line 232

define_getter :flag

#index_loaded? ⇒ Boolean

Returns:

  • (Boolean) 


145
146
147
148
149
 
# File 'lib/hts/bam.rb', line 145

def index_loaded?
  check_closed

  @native.index_loaded?
end

#insert_size ⇒ Array Also known as: isize

Get insert_size array

Returns:

  • (Array) —

    the insert_size array



239
 
# File 'lib/hts/bam.rb', line 239

define_getter :insert_size

#load_index(index_name = nil) ⇒ Object

Raises:



131
132
133
134
135
 
# File 'lib/hts/bam.rb', line 131

def load_index(index_name = nil)
  return self if try_load_index(index_name)

  raise MissingIndexError, "Failed to load index #{index_name || "for #{@file_name}"}"
end

#mapq ⇒ Array

Get mapq array

Returns:

  • (Array) —

    the mapq array



235
 
# File 'lib/hts/bam.rb', line 235

define_getter :mapq

#mate_chrom ⇒ Array

Get mate_chrom array

Returns:

  • (Array) —

    the mate_chrom array



237
 
# File 'lib/hts/bam.rb', line 237

define_getter :mate_chrom

#mate_pos ⇒ Array Also known as: mpos

Get mate_pos array

Returns:

  • (Array) —

    the mate_pos array



238
 
# File 'lib/hts/bam.rb', line 238

define_getter :mate_pos

#pileup(region = nil, beg = nil, end_: nil, maxcnt: nil, &block) ⇒ Object

Pileup iterator over this file. Optional region can be specified. When a block is given, uses RAII-style and ensures the iterator is closed at block end. Without a block, returns an Enumerator over a live Pileup instance; caller should close when done.

Parameters:

  • region (String, nil) (defaults to: nil) —

    region string like "chr1:100-200"

  • beg (Integer, nil) (defaults to: nil)
  • end_ (Integer, nil) (defaults to: nil)
  • maxcnt (Integer, nil) (defaults to: nil) —

    cap on depth per position



362
363
364
365
366
367
368
369
370
371
372
373
 
# File 'lib/hts/bam.rb', line 362

def pileup(region = nil, beg = nil, end_: nil, maxcnt: nil, &block)
  check_closed
  if block_given?
    Pileup.open(self, region:, beg:, end_: end_, maxcnt: maxcnt) do |piter|
      piter.each(&block)
    end
    self
  else
    piter = Pileup.new(self, region:, beg:, end_: end_, maxcnt: maxcnt)
    piter.to_enum(:each)
  end
end

#pos ⇒ Array

Get pos array

Returns:

  • (Array) —

    the pos array



234
 
# File 'lib/hts/bam.rb', line 234

define_getter :pos

#qname ⇒ Array

Get qname array

Returns:

  • (Array) —

    the qname array



231
 
# File 'lib/hts/bam.rb', line 231

define_getter :qname

#qual ⇒ Array

Get qual array

Returns:

  • (Array) —

    the qual array



241
 
# File 'lib/hts/bam.rb', line 241

define_getter :qual

#query(region, beg = nil, end_ = nil, copy: false, &block) ⇒ Object

Iterate records in a genomic region or multiple regions. See #each for copy semantics. When copy: false, the yielded Record is reused and should not be stored.

Examples:

Single region query

bam.query("chr1:100-200") { |r| puts r.qname }
bam.query("chr1", 100, 200) { |r| puts r.qname }

Multi-region query

bam.query(["chr1:100-200", "chr2:500-600"]) { |r| puts r.qname }

Parameters:

  • region (String, Array<String>) —

    Region specification(s)

    • Single region: "chr1:100-200" or "chr1" with beg/end parameters
    • Multiple regions: ["chr1:100-200", "chr2:500-600", ...]
  • beg (Integer, nil) (defaults to: nil) —

    Start position (used with single string region)

  • end_ (Integer, nil) (defaults to: nil) —

    End position (used with single string region)

  • copy (Boolean) (defaults to: false) —

    Whether to deep-copy records (see #each)



331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
 
# File 'lib/hts/bam.rb', line 331

def query(region, beg = nil, end_ = nil, copy: false, &block)
  check_closed
  raise "Index file is required to call the query method." unless ensure_index_loaded

  case region
  when Array
    raise ArgumentError, "beg and end_ cannot be used with array of regions" if beg || end_

    query_regions(region, copy:, &block)
  when String
    if beg && end_
      tid = header.get_tid(region)
      queryi(tid, beg, end_, copy:, &block)
    elsif beg.nil? && end_.nil?
      querys(region, copy:, &block)
    else
      raise ArgumentError, "beg and end_ must be specified together"
    end
  else
    raise ArgumentError, "region must be String or Array"
  end
end

#rewind ⇒ Object 



184
185
186
187
188
189
190
191
 
# File 'lib/hts/bam.rb', line 184

def rewind
  raise "Cannot rewind: no start position" unless @start_position

  result = seek(@start_position)
  raise "Failed to rewind: #{result}" if result.negative?

  tell
end

#seek(offset) ⇒ Object 



181
 
# File 'lib/hts/bam.rb', line 181

def seek(offset) = @native.seek(offset)

#seq ⇒ Array

Get seq array

Returns:

  • (Array) —

    the seq array



240
 
# File 'lib/hts/bam.rb', line 240

define_getter :seq

#set_threads(n = nil) ⇒ Object

Raises:

  • (TypeError) 


168
169
170
171
172
173
174
175
176
177
178
179
 
# File 'lib/hts/bam.rb', line 168

def set_threads(n = nil)
  if n.nil?
    require "etc"
    n = [Etc.nprocessors - 1, 1].max
  end
  raise TypeError unless n.is_a?(Integer)
  raise ArgumentError, "Number of threads must be positive" if n < 1
  raise "Failed to set number of threads: #{n}" if @native.set_threads(n).negative?

  @nthreads = n
  self
end

#tell ⇒ Object 



182
 
# File 'lib/hts/bam.rb', line 182

def tell = @native.tell

#to_a ⇒ Object

Materialize independent, owning records. Enumerable#to_a is unsafe for the default reused-record iterator because every array element would otherwise refer to the same native bam1_t buffer.



266
267
268
 
# File 'lib/hts/bam.rb', line 266

def to_a
  collect_records
end

#try_load_index(index_name = nil) ⇒ Object 



137
138
139
140
141
142
143
 
# File 'lib/hts/bam.rb', line 137

def try_load_index(index_name = nil)
  check_closed

  @index_name = index_name
  @index_load_attempted = true
  @native.load_index(index_name)
end

#write(record) ⇒ Object 



213
214
215
216
217
218
219
220
 
# File 'lib/hts/bam.rb', line 213

def write(record)
  check_closed

  r = @native.write(header.__send__(:native_handle), record.__send__(:native_handle))
  raise "Failed to write record" if r < 0

  nil
end

#write_header(header) ⇒ Object 



201
202
203
204
205
206
207
 
# File 'lib/hts/bam.rb', line 201

def write_header(header)
  check_closed

  @header = header.dup
  @native.write_header(header.__send__(:native_handle))
  nil
end