Class: HTS::Bam::Pileup

Inherits:
Object
  • Object
show all
Includes:
Enumerable
Defined in:
lib/hts/bam/pileup.rb

Overview

High-level pileup iterator for a single SAM/BAM/CRAM

Defined Under Namespace

Classes: PileupColumn, PileupRecord

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(bam, region: nil, beg: nil, end_: nil, maxcnt: nil) ⇒ Pileup

Create a Pileup iterator

Parameters:

  • bam (HTS::Bam)
  • region (String, nil) (defaults to: nil) —

    Optional region string (requires index)

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

    Optional begin when using tid/beg/end form

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

    Optional end when using tid/beg/end form

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

    Max per-position depth (capped)



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

def initialize(bam, region: nil, beg: nil, end_: nil, maxcnt: nil)
  @bam    = bam
  @header = bam.header
  @itr    = nil
  @cb     = nil
  @plp    = nil

  # Optional region iterator
  if region && beg.nil? && end_.nil?
    raise "Index file is required to use region pileup." unless bam.index_loaded?

    @itr = HTS::LibHTS.sam_itr_querys(bam.instance_variable_get(:@idx), @header.struct, region)
    raise "Failed to query region: #{region}" if @itr.null?
  elsif region && beg && end_
    raise "Index file is required to use region pileup." unless bam.index_loaded?

    tid = @header.get_tid(region)
    @itr = HTS::LibHTS.sam_itr_queryi(bam.instance_variable_get(:@idx), tid, beg, end_)
    raise "Failed to query region: #{region} #{beg} #{end_}" if @itr.null?
  elsif beg || end_
    raise ArgumentError, "beg and end_ must be specified together"
  end

  # Build the auto callback for bam_plp_init (micro-optimized)
  # - Hoist ivar/constant lookups out of the callback to reduce per-call overhead.
  # - Specialize callbacks to avoid branching in the hot path.
  hts_fp     = @bam.instance_variable_get(:@hts_file)
  hdr_struct = @header.struct
  itr_local  = @itr

  @cb = if itr_local && !itr_local.null?
          FFI::Function.new(:int, %i[pointer pointer]) do |_data, b|
            # HTSlib contract: sam_itr_next returns >= 0 on success, -1 on EOF, < -1 on error.
            r = HTS::LibHTS.sam_itr_next(hts_fp, itr_local, b)
            if r >= 0
              0
            else
              (r == -1 ? -1 : -2)
            end
          end
        else
          FFI::Function.new(:int, %i[pointer pointer]) do |_data, b|
            # HTSlib contract: sam_read1 returns >= 0 on success, -1 on EOF/error.
            r = HTS::LibHTS.sam_read1(hts_fp, hdr_struct, b)
            r == -1 ? -1 : 0
          end
        end

  @plp = HTS::LibHTS.bam_plp_init(@cb, nil)
  raise "bam_plp_init failed" if @plp.null?

  HTS::LibHTS.bam_plp_set_maxcnt(@plp, maxcnt) if maxcnt
end

Class Method Details

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

Usage:

HTS::Bam::Pileup.open(bam, region: "chr1:1-100") do |pl|
  pl.each { |col| ... }
end



13
14
15
16
17
18
19
20
21
22
23
 
# File 'lib/hts/bam/pileup.rb', line 13

def self.open(*args, **kw)
  pu = new(*args, **kw)
  return pu unless block_given?

  begin
    yield pu
  ensure
    pu.close
  end
  pu
end

Instance Method Details

#close ⇒ Object 



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

def close
  if @plp && !@plp.null?
    HTS::LibHTS.bam_plp_destroy(@plp)
    @plp = FFI::Pointer::NULL
  end
  if @itr && !@itr.null?
    HTS::LibHTS.hts_itr_destroy(@itr)
    @itr = FFI::Pointer::NULL
  end
  # Keep @cb referenced by instance to avoid GC during iteration.
  @cb
end

#each ⇒ Object 



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

def each
  return to_enum(__method__) unless block_given?

  tid_ptr = FFI::MemoryPointer.new(:int)
  pos_ptr = FFI::MemoryPointer.new(:long_long) # hts_pos_t
  n_ptr   = FFI::MemoryPointer.new(:int)

  # Micro-optimizations:
  # - Compute constant struct size once
  # - Hoist header reference outside the loop
  plp1_size    = HTS::LibHTS::BamPileup1.size
  header_local = @header

  loop do
    base_ptr = HTS::LibHTS.bam_plp64_auto(@plp, tid_ptr, pos_ptr, n_ptr)

    # When base_ptr is NULL, check n to distinguish EOF (n == 0) from error (n < 0)
    if base_ptr.null?
      n = n_ptr.read_int
      raise "HTSlib pileup error (bam_plp64_auto)" if n < 0

      break
    end

    tid = tid_ptr.read_int
    pos = pos_ptr.read_long_long
    n   = n_ptr.read_int

    # Construct alignment entries with minimal allocations
    if n.zero?
      alignments = []
    else
      alignments = Array.new(n)
      i = 0
      while i < n
        e_ptr = base_ptr + (i * plp1_size)
        entry = HTS::LibHTS::BamPileup1.new(e_ptr)
        alignments[i] = PileupRecord.new(entry, header_local)
        i += 1
      end
    end

    yield PileupColumn.new(tid: tid, pos: pos, alignments: alignments)
  end

  self
end

#reset ⇒ Object 



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

def reset
  HTS::LibHTS.bam_plp_reset(@plp) if @plp && !@plp.null?
end