Class: HTS::Bam::Mpileup

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

Overview

High-level mpileup iterator over multiple BAM/CRAM inputs

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(inputs, region: nil, beg: nil, end_: nil, maxcnt: nil, overlaps: false) ⇒ Mpileup

Normalize inputs to HTS::Bam instances Accepts array of HTS::Bam or filenames (String)

Raises:

  • (ArgumentError) 


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

def initialize(inputs, region: nil, beg: nil, end_: nil, maxcnt: nil, overlaps: false)
  raise ArgumentError, "inputs must be non-empty" if inputs.nil? || inputs.empty?

  @owned_bams = [] # Bams we opened here; will be closed on close
  @bams = inputs.map do |x|
    case x
    when HTS::Bam
      x
    when String
      b = HTS::Bam.open(x)
      @owned_bams << b
      b
    else
      raise ArgumentError, "Unsupported input type: #{x.class}"
    end
  end

  n = @bams.length
  @iters       = []
  @data_blocks = [] # per-input packed pointers kept alive

  # Prepare optional region iterators for each input
  @bams.each_with_index do |bam, i|
    itr = nil
    if region && beg.nil? && end_.nil?
      raise "Index required for region mpileup" unless bam.index_loaded?

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

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

  # Build per-input packed pointer blocks so C passes them back to the callback.
  # Layout per input: [0] hts_fp (htsFile*), [1] hdr_struct (bam_hdr_t*), [2] itr (hts_itr_t* or NULL)
  ptr_size = FFI.type_size(:pointer)
  data_array = FFI::MemoryPointer.new(:pointer, n)
  @bams.each_with_index do |bam, i|
    hts_fp     = bam.instance_variable_get(:@hts_file)
    hdr_struct = bam.header.struct
    itr        = @iters[i]
    block = FFI::MemoryPointer.new(:pointer, 3)
    block.put_pointer(0 * ptr_size, hts_fp)
    block.put_pointer(1 * ptr_size, hdr_struct)
    block.put_pointer(2 * ptr_size, itr && !itr.null? ? itr : FFI::Pointer::NULL)
    @data_blocks << block
    data_array.put_pointer(i * ptr_size, block)
  end
  # Keep the array of per-input blocks alive while the C side holds on to them
  @data_array = data_array

  @cb = FFI::Function.new(:int, %i[pointer pointer]) do |data, b|
    # Unpack pointers from the per-input block
    hts_fp     = data.get_pointer(0 * ptr_size)
    hdr_struct = data.get_pointer(1 * ptr_size)
    itr        = data.get_pointer(2 * ptr_size)
    if itr && !itr.null?
      r = HTS::LibHTS.sam_itr_next(hts_fp, itr, b)
      if r >= 0
        0
      else
        (r == -1 ? -1 : -2)
      end
    else
      r = HTS::LibHTS.sam_read1(hts_fp, hdr_struct, b)
      r == -1 ? -1 : 0
    end
  end

  @iter = HTS::LibHTS.bam_mplp_init(n, @cb, @data_array)
  raise "bam_mplp_init failed" if @iter.null?

  HTS::LibHTS.bam_mplp_set_maxcnt(@iter, maxcnt) if maxcnt
  return unless overlaps

  rc = HTS::LibHTS.bam_mplp_init_overlaps(@iter)
  raise "bam_mplp_init_overlaps failed" if rc < 0
end

Class Method Details

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

Usage:

HTS::Bam::Mpileup.open([bam1, bam2], region: "chr1:1-100") do |mpl|
  mpl.each { |cols| ... }
end



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

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

  begin
    yield m
  ensure
    m.close
  end
  m
end

Instance Method Details

#close ⇒ Object 



160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
 
# File 'lib/hts/bam/mpileup.rb', line 160

def close
  if @iter && !@iter.null?
    HTS::LibHTS.bam_mplp_destroy(@iter)
    @iter = FFI::Pointer::NULL
  end
  @iters.each do |itr|
    HTS::LibHTS.hts_itr_destroy(itr) if itr && !itr.null?
  end
  @iters.clear
  # Keep references to callback and data blocks to prevent GC
  @_keepalive = [@cb, @data_array, *@data_blocks]
  # Close owned bams opened by this object
  @owned_bams.each do |b|
    b.close
  rescue StandardError
  end
  @owned_bams.clear
end

#each ⇒ Object

Yields an array of Pileup::PileupColumn (one per input) for each position



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

def each
  return to_enum(__method__) unless block_given?

  n = @bams.length
  tid_ptr = FFI::MemoryPointer.new(:int)
  pos_ptr = FFI::MemoryPointer.new(:long_long)
  n_ptr   = FFI::MemoryPointer.new(:int, n)
  plp_ptr = FFI::MemoryPointer.new(:pointer, n)
  plp1_size = HTS::LibHTS::BamPileup1.size
  headers   = @bams.map(&:header)

  while HTS::LibHTS.bam_mplp64_auto(@iter, tid_ptr, pos_ptr, n_ptr, plp_ptr) > 0
    tid = tid_ptr.read_int
    pos = pos_ptr.read_long_long

    counts = n_ptr.read_array_of_int(n)
    plp_arr = plp_ptr.read_array_of_pointer(n)

    cols = Array.new(n)
    i = 0
    while i < n
      c = counts[i]
      if c <= 0 || plp_arr[i].null?
        cols[i] = HTS::Bam::Pileup::PileupColumn.new(tid: tid, pos: pos, alignments: [])
      else
        base_ptr = plp_arr[i]
        aligns = Array.new(c)
        j = 0
        while j < c
          e_ptr = base_ptr + (j * plp1_size)
          entry = HTS::LibHTS::BamPileup1.new(e_ptr)
          aligns[j] = HTS::Bam::Pileup::PileupRecord.new(entry, headers[i])
          j += 1
        end
        cols[i] = HTS::Bam::Pileup::PileupColumn.new(tid: tid, pos: pos, alignments: aligns)
      end
      i += 1
    end

    yield cols
  end

  self
end