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

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)



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

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
  @cb = FFI::Function.new(:int, %i[pointer pointer]) do |_data, b|
    if @itr && !@itr.null?
      slen = HTS::LibHTS.sam_itr_next(@bam.instance_variable_get(:@hts_file), @itr, b)
      if slen > 0
        0
      elsif slen == -1
        -1
      else
        -2
      end
    else
      r = HTS::LibHTS.sam_read1(@bam.instance_variable_get(:@hts_file), @header.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

Instance Method Details

#close ⇒ Object 



156
157
158
159
160
161
162
163
164
165
166
167
 
# File 'lib/hts/bam/pileup.rb', line 156

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 



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

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)

  begin
    while (base_ptr = HTS::LibHTS.bam_plp64_auto(@plp, tid_ptr, pos_ptr, n_ptr)) && !base_ptr.null?
      tid = tid_ptr.read_int
      pos = pos_ptr.read_long_long
      n   = n_ptr.read_int

      # Construct alignment entries
      alignments = if n.zero?
                     []
                   else
                     size = HTS::LibHTS::BamPileup1.size
                     n.times.map do |i|
                       e_ptr = base_ptr + (i * size)
                       entry = HTS::LibHTS::BamPileup1.new(e_ptr)
                       PileupRecord.new(entry, @header)
                     end
                   end

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

  self
end

#reset ⇒ Object 



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

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