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

Methods inherited from Hts

#closed?, #fai=, #file_format, #file_format_version, #rewind, #seek, #set_threads, #struct, #tell, #to_ptr

Constructor Details

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

Returns a new instance of Bam.

Raises:

  • (Errno::ENOENT) 


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

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
  @hts_file   = LibHTS.hts_open(@file_name, mode)

  raise Errno::ENOENT, "Failed to open #{@file_name}" if @hts_file.null?

  # 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 = LibHTS.hts_set_fai_filename(@hts_file, 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(@hts_file)
  build_index(index) if build_index
  @idx = load_index(index)
  @start_position = tell
end

Instance Attribute Details

#file_name ⇒ Object (readonly)

Returns the value of attribute file_name.



20
21
22
 
# File 'lib/hts/bam.rb', line 20

def file_name
  @file_name
end

#header ⇒ Object

Returns the value of attribute header.



20
21
22
 
# File 'lib/hts/bam.rb', line 20

def header
  @header
end

#index_name ⇒ Object (readonly)

Returns the value of attribute index_name.



20
21
22
 
# File 'lib/hts/bam.rb', line 20

def index_name
  @index_name
end

#mode ⇒ Object (readonly)

Returns the value of attribute mode.



20
21
22
 
# File 'lib/hts/bam.rb', line 20

def mode
  @mode
end

#nthreads ⇒ Object (readonly)

Returns the value of attribute nthreads.



20
21
22
 
# File 'lib/hts/bam.rb', line 20

def nthreads
  @nthreads
end

Class Method Details

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



34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
 
# File 'lib/hts/bam.rb', line 34

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 LibHTS.sam_index_build3(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

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



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

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 



143
144
145
 
# File 'lib/hts/bam.rb', line 143

def <<(record)
  write(record)
end

#aux(tag) ⇒ Object

FIXME: experimental



167
168
169
170
171
172
173
174
 
# File 'lib/hts/bam.rb', line 167

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 



96
97
98
99
100
101
 
# File 'lib/hts/bam.rb', line 96

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)
  self # for method chaining
end

#chrom ⇒ Array

Get chrom array

Returns:

  • (Array) —

    the chrom array



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

define_getter :chrom

#cigar ⇒ Array

Get cigar array

Returns:

  • (Array) —

    the cigar array



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

define_getter :cigar

#close ⇒ Object 



119
120
121
122
123
 
# File 'lib/hts/bam.rb', line 119

def close
  LibHTS.hts_idx_destroy(@idx) if @idx && !@idx.null?
  @idx = nil
  super
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.



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

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

#each_aux(tag) ⇒ Object

FIXME: experimental



195
196
197
198
199
200
201
202
203
204
 
# File 'lib/hts/bam.rb', line 195

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



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

define_iterator :chrom

#each_cigar ⇒ Object

Get cigar iterator



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

define_iterator :cigar

#each_flag ⇒ Object

Get flag iterator



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

define_iterator :flag

#each_insert_size ⇒ Object Also known as: each_isize

Get insert_size iterator



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

define_iterator :insert_size

#each_mapq ⇒ Object

Get mapq iterator



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

define_iterator :mapq

#each_mate_chrom ⇒ Object

Get mate_chrom iterator



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

define_iterator :mate_chrom

#each_mate_pos ⇒ Object Also known as: each_mpos

Get mate_pos iterator



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

define_iterator :mate_pos

#each_pos ⇒ Object

Get pos iterator



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

define_iterator :pos

#each_qname ⇒ Object

Get qname iterator



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

define_iterator :qname

#each_qual ⇒ Object

Get qual iterator



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

define_iterator :qual

#each_seq ⇒ Object

Get seq iterator



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

define_iterator :seq

#flag ⇒ Array

Get flag array

Returns:

  • (Array) —

    the flag array



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

define_getter :flag

#index_loaded? ⇒ Boolean

Returns:

  • (Boolean) 


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

def index_loaded?
  check_closed

  !@idx.null?
end

#insert_size ⇒ Array Also known as: isize

Get insert_size array

Returns:

  • (Array) —

    the insert_size array



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

define_getter :insert_size

#load_index(index_name = nil) ⇒ Object 



103
104
105
106
107
108
109
110
111
 
# File 'lib/hts/bam.rb', line 103

def load_index(index_name = nil)
  check_closed

  if index_name
    LibHTS.sam_index_load2(@hts_file, @file_name, index_name)
  else
    LibHTS.sam_index_load3(@hts_file, @file_name, nil, 2) # should be 3 ? (copy remote file to local?)
  end
end

#mapq ⇒ Array

Get mapq array

Returns:

  • (Array) —

    the mapq array



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

define_getter :mapq

#mate_chrom ⇒ Array

Get mate_chrom array

Returns:

  • (Array) —

    the mate_chrom array



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

define_getter :mate_chrom

#mate_pos ⇒ Array Also known as: mpos

Get mate_pos array

Returns:

  • (Array) —

    the mate_pos array



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

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



268
269
270
271
272
273
274
275
276
277
278
279
 
# File 'lib/hts/bam.rb', line 268

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



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

define_getter :pos

#qname ⇒ Array

Get qname array

Returns:

  • (Array) —

    the qname array



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

define_getter :qname

#qual ⇒ Array

Get qual array

Returns:

  • (Array) —

    the qual array



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

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)



237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
 
# File 'lib/hts/bam.rb', line 237

def query(region, beg = nil, end_ = nil, copy: false, &block)
  check_closed
  raise "Index file is required to call the query method." unless 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

#seq ⇒ Array

Get seq array

Returns:

  • (Array) —

    the seq array



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

define_getter :seq

#write(record) ⇒ Object 



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

def write(record)
  check_closed

  r = LibHTS.sam_write1(@hts_file, header, record)
  raise "Failed to write record" if r < 0
end

#write_header(header) ⇒ Object 



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

def write_header(header)
  check_closed

  @header = header.dup
  LibHTS.sam_hdr_write(@hts_file, header)
end