Class: Bio::DB::Sam

Inherits:
Object
  • Object
show all
Defined in:
lib/bio/db/sam.rb

Instance Attribute Summary collapse

Class Method Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(optsa = {}) ⇒ Sam

Returns a new instance of Sam.



18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 
# File 'lib/bio/db/sam.rb', line 18

def initialize(optsa={})
  opts =  { :fasta => nil,  :bam => nil,:tam => nil, :compressed => true, :write => false }.merge!(optsa)



  @fasta_path = opts[:fasta]
  @compressed = opts[:compressed]
  @write = opts[:write]
  bam = opts[:bam]
  tam = opts[:tam]

  if bam == nil && tam == nil && @fasta_path == nil then
    raise SAMException.new(), "No alignment or reference file"
  elsif bam != nil && tam != nil then
    raise SAMException.new(), "Alignment has to be in either text or binary format, not both"
  elsif bam != nil then
    @binary = true
    @sam = bam     
  elsif tam != nil then
    @sam = tam     
    @binary = false

  end
  @fasta_file = nil
  @sam_file   = nil

  ObjectSpace.define_finalizer(self,  self.class.method(:finalize).to_proc)
end

Instance Attribute Details

#sam_file ⇒ Object (readonly)

Returns the value of attribute sam_file.



16
17
18
 
# File 'lib/bio/db/sam.rb', line 16

def sam_file
  @sam_file
end

Class Method Details

.finalize(id) ⇒ Object 



87
88
89
90
 
# File 'lib/bio/db/sam.rb', line 87

def Sam.finalize(id)
  id.close()
  puts "Finalizing #{id}  at #{Time.new}"       
end

Instance Method Details

#average_coverage(chromosome, qstart, len) ⇒ Object 



117
118
119
120
121
122
123
124
125
126
127
128
129
 
# File 'lib/bio/db/sam.rb', line 117

def average_coverage(chromosome, qstart, len)

  #reference = fetch_reference(chromosome, qstart,len)
  # len = reference.length if len > reference.length


  coverages = chromosome_coverage(chromosome, qstart, len)
  total = 0
  len.times{ |i| total= total + coverages[i] }
  avg_cov = total.to_f / len
  #LibC.free reference
  avg_cov
end

#chromosome_coverage(chromosome, qstart, len) ⇒ Object 



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
159
160
161
 
# File 'lib/bio/db/sam.rb', line 131

def chromosome_coverage(chromosome, qstart, len)
  #  reference = fetch_reference(chromosome, qstart,len)
  #  len = reference.length if len > reference.length
  #p qend.to_s + "-" + qstart.to_s + "framesize " + (qend - qstart).to_s
  coverages = Array.new(len, 0)

  chr_cov_proc = Proc.new do |alignment|
    #last = qstart + len
    #first = qstart
    #last = alignment.calend if last > alignment.calend
    #first = alignment.pos if first < alignment.pos
    # p first
    last = alignment.calend - qstart
    first = alignment.pos - qstart
    if last < first
      tmp = last
      last = first
      first = last
    end

    #  STDERR.puts "#{first} #{last}\n"
    first.upto(last-1) { |i|

      coverages[i-1] = 1 + coverages[i-1]  if i-1 < len && i > 0
    }
  end

  fetch_with_function(chromosome, qstart, qstart+len,  chr_cov_proc)
  #p coverages
  coverages
end

#close ⇒ Object 



79
80
81
82
83
84
85
 
# File 'lib/bio/db/sam.rb', line 79

def close()
  Bio::DB::SAM::Tools.fai_destroy(@fasta_index) unless @fasta_index.nil? || @fasta_index.null?
  Bio::DB::SAM::Tools.bam_index_destroy(@sam_index) unless @sam_index.nil? || @sam_index.null?
  Bio::DB::SAM::Tools.samclose(@sam_file) unless @sam_file.nil? 
  @sam_file = nil
  @fasta_index = nil
end

#fetch(chromosome, qstart, qend) ⇒ Object 



178
179
180
181
182
183
184
185
186
 
# File 'lib/bio/db/sam.rb', line 178

def fetch(chromosome, qstart, qend)
  als = Array.new
  fetchAlignment = Proc.new do |alignment|
    als.push(alignment.clone)   
    0  
  end
  fetch_with_function(chromosome, qstart, qend, fetchAlignment)
  als
end

#fetch_reference(chromosome, qstart, qend) ⇒ Object

Raises:



163
164
165
166
167
168
169
170
171
 
# File 'lib/bio/db/sam.rb', line 163

def fetch_reference(chromosome, qstart,qend)
  load_reference if @fasta_index.nil? || @fasta_index.null? 
  query = query_string(chromosome, qstart,qend)
  len = FFI::MemoryPointer.new :int
  reference = Bio::DB::SAM::Tools.fai_fetch(@fasta_index, query, len)
  raise SAMException.new(), "Unable to get sequence for reference: "+query if reference.nil?

  reference
end

#fetch_with_function(chromosome, qstart, qend, function) ⇒ Object 



188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
 
# File 'lib/bio/db/sam.rb', line 188

def fetch_with_function(chromosome, qstart, qend, function)
  load_index if @sam_index.nil? || @sam_index.null?
  chr = FFI::MemoryPointer.new :int
  beg = FFI::MemoryPointer.new :int
  last = FFI::MemoryPointer.new :int
  query = query_string(chromosome, qstart,qend)
  qpointer = FFI::MemoryPointer.from_string(query)
  header = @sam_file[:header]
  Bio::DB::SAM::Tools.bam_parse_region(header,qpointer, chr, beg, last) 
  #raise SAMException.new(), "invalid query: " + query  if(chr.read_int < 0)
  count = 0;

  fetchAlignment = Proc.new do |bam_alignment, data|
    alignment =  Alignment.new
    alignment.set(bam_alignment, header)
    function.call(alignment)
    count = count + 1
    0  
  end
  Bio::DB::SAM::Tools.bam_fetch(@sam_file[:x][:bam], @sam_index,chr.read_int,beg.read_int, last.read_int, nil, fetchAlignment)
  #LibC.free chr
  #LibC.free beg
  #LibC.free last
  #LibC.free qpointer
  count
end

#load_index ⇒ Object

Raises:



92
93
94
95
96
97
98
99
100
101
 
# File 'lib/bio/db/sam.rb', line 92

def load_index()
  raise SAMException.new(), "Indexes are only supported by BAM files, please use samtools to convert your SAM file" unless @binary
  @sam_index = Bio::DB::SAM::Tools.bam_index_load(@sam)
  if @sam_index.null? then
    p "Generating index for: " + @sam
    Bio::DB::SAM::Tools.bam_index_build(@sam)
    @sam_index = Bio::DB::SAM::Tools.bam_index_load(@sam)
    raise SAMException.new(), "Unable to generate bam index for: " + @sam if @sam_index.nil? || @sam_index.null?
  end
end

#load_reference ⇒ Object

Raises:



103
104
105
106
107
108
109
110
111
112
113
114
115
 
# File 'lib/bio/db/sam.rb', line 103

def load_reference()
  raise SAMException.new(), "No path for the refernce fasta file. " if @fasta_path.nil?

  @fasta_index = Bio::DB::SAM::Tools.fai_load(@fasta_path)

  if @fasta_index.null? then
    p "Generating index for: " + @fasta_path
    Bio::DB::SAM::Tools.fai_build(@fasta_path)
    @fasta_index =  Bio::DB::SAM::Tools.fai_load(@fasta_path)
    raise SAMException.new(), "Unable to generate fasta index for: " + @fasta_path if @fasta_index.nil? ||  @fasta_index.null?
  end

end

#open ⇒ Object

Raises:



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
 
# File 'lib/bio/db/sam.rb', line 47

def open()

  raise SAMException.new(), "Writing not supported yet" if @write
  raise SAMException.new(), "No SAM file specified" unless @sam 

  opts = @write ? "w" : "r"
  if @binary then  
    opts += "b" 
    if @write then
      unless @compressed then 
        opts += "u"
      end
    end
  end
  valid = ["r", "w", "wh", "rb", "wb" , "wbu"]
  unless valid.include?(opts) then
    raise SAMException.new(), "Invalid options for samopen: " + opts 
  end

  samFile = Bio::DB::SAM::Tools.samopen(@sam, opts, nil)
  if samFile.null? then
    @sam_file = nil
    raise SAMException.new(), "File not opened:  " + @sam
  end
  @sam_file = Bio::DB::SAM::Tools::SamfileT.new(samFile)

end

#query_string(chromosome, qstart, qend) ⇒ Object 



173
174
175
176
 
# File 'lib/bio/db/sam.rb', line 173

def query_string(chromosome, qstart,qend)
  query = chromosome + ":" + qstart.to_s + "-" + qend.to_s 
  query
end

#to_s ⇒ Object 



75
76
77
 
# File 'lib/bio/db/sam.rb', line 75

def to_s()
  (@binary ? "Binary" : "Text") + " file: " + @sam + " with fasta: " + @fasta_path
end