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, Mpileup, Pileup, Record

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.



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

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

  return if @mode[0] == "w"

  @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.



34
35
36
 
# File 'lib/hts/bam.rb', line 34

def file_name
  @file_name
end

#header ⇒ Object

Returns the value of attribute header.



34
35
36
 
# File 'lib/hts/bam.rb', line 34

def header
  @header
end

#index_name ⇒ Object (readonly)

Returns the value of attribute index_name.



34
35
36
 
# File 'lib/hts/bam.rb', line 34

def index_name
  @index_name
end

#mode ⇒ Object (readonly)

Returns the value of attribute mode.



34
35
36
 
# File 'lib/hts/bam.rb', line 34

def mode
  @mode
end

#nthreads ⇒ Object (readonly)

Returns the value of attribute nthreads.



34
35
36
 
# File 'lib/hts/bam.rb', line 34

def nthreads
  @nthreads
end

Class Method Details

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



48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
 
# File 'lib/hts/bam.rb', line 48

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.



21
22
23
24
25
26
27
28
29
30
31
32
 
# File 'lib/hts/bam.rb', line 21

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 



36
37
38
39
40
41
42
43
44
45
46
 
# File 'lib/hts/bam.rb', line 36

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

  begin
    yield file
  ensure
    file.close
  end
  file
end

Instance Method Details

#<<(record) ⇒ Object 



193
194
195
 
# File 'lib/hts/bam.rb', line 193

def <<(record)
  write(record)
end

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

FIXME: experimental



217
218
219
220
221
222
223
224
 
# File 'lib/hts/bam.rb', line 217

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 



113
114
115
116
117
118
119
120
 
# File 'lib/hts/bam.rb', line 113

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



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

define_getter :chrom

#cigar ⇒ Array

Get cigar array

Returns:

  • (Array) —

    the cigar array



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

define_getter :cigar

#close ⇒ Object 



136
137
138
 
# File 'lib/hts/bam.rb', line 136

def close
  @native&.close
end

#closed? ⇒ Boolean

Returns:

  • (Boolean) 


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

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.



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

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.


270
271
272
273
274
275
276
 
# File 'lib/hts/bam.rb', line 270

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

#each_aux(tag) ⇒ Object

FIXME: experimental



252
253
254
255
256
257
258
259
260
261
 
# File 'lib/hts/bam.rb', line 252

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



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

define_iterator :chrom

#each_cigar ⇒ Object

Get cigar iterator



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

define_iterator :cigar

#each_flag ⇒ Object

Get flag iterator



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

define_iterator :flag

#each_insert_size ⇒ Object Also known as: each_isize

Get insert_size iterator



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

define_iterator :insert_size

#each_mapq ⇒ Object

Get mapq iterator



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

define_iterator :mapq

#each_mate_chrom ⇒ Object

Get mate_chrom iterator



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

define_iterator :mate_chrom

#each_mate_pos ⇒ Object Also known as: each_mpos

Get mate_pos iterator



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

define_iterator :mate_pos

#each_pos ⇒ Object

Get pos iterator



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

define_iterator :pos

#each_qname ⇒ Object

Get qname iterator



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

define_iterator :qname

#each_qual ⇒ Object

Get qual iterator



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

define_iterator :qual

#each_seq ⇒ Object

Get seq iterator



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

define_iterator :seq

#file_format ⇒ Object 



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

def file_format = @native.file_format

#file_format_version ⇒ Object 



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

def file_format_version = @native.file_format_version

#flag ⇒ Array

Get flag array

Returns:

  • (Array) —

    the flag array



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

define_getter :flag

#index_loaded? ⇒ Boolean

Returns:

  • (Boolean) 


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

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



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

define_getter :insert_size

#load_index(index_name = nil) ⇒ Object 



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

def load_index(index_name = nil)
  check_closed

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

#mapq ⇒ Array

Get mapq array

Returns:

  • (Array) —

    the mapq array



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

define_getter :mapq

#mate_chrom ⇒ Array

Get mate_chrom array

Returns:

  • (Array) —

    the mate_chrom array



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

define_getter :mate_chrom

#mate_pos ⇒ Array Also known as: mpos

Get mate_pos array

Returns:

  • (Array) —

    the mate_pos array



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

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



325
326
327
328
329
330
331
332
333
334
335
336
 
# File 'lib/hts/bam.rb', line 325

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



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

define_getter :pos

#qname ⇒ Array

Get qname array

Returns:

  • (Array) —

    the qname array



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

define_getter :qname

#qual ⇒ Array

Get qual array

Returns:

  • (Array) —

    the qual array



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

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)



294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
 
# File 'lib/hts/bam.rb', line 294

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 



160
161
162
163
164
165
166
167
 
# File 'lib/hts/bam.rb', line 160

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 



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

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

#seq ⇒ Array

Get seq array

Returns:

  • (Array) —

    the seq array



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

define_getter :seq

#set_threads(n = nil) ⇒ Object

Raises:

  • (TypeError) 


144
145
146
147
148
149
150
151
152
153
154
155
 
# File 'lib/hts/bam.rb', line 144

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 



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

def tell = @native.tell

#write(record) ⇒ Object 



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

def write(record)
  check_closed

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

#write_header(header) ⇒ Object 



175
176
177
178
179
180
 
# File 'lib/hts/bam.rb', line 175

def write_header(header)
  check_closed

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