Permalink
Browse files

-Added the IDs to test

-Capacity to calculate the average coverage of certain region without unnecesary allocation.
  • Loading branch information...
1 parent 91a67ac commit fd0cd996e5dc477c3e6283f12e5c4abb3ae933be @homonecloco committed Aug 23, 2010
Showing with 70 additions and 7 deletions.
  1. +7 −6 Rakefile
  2. +27 −1 lib/bio/db/sam.rb
  3. +35 −0 test/basictest.rb
  4. +1 −0 test/samples/small/ids2.txt
View
@@ -1,10 +1,11 @@
-require 'echoe'
+#require 'echoe'
-Echoe.new("samtools-ruby", version=0.1) do |p|
- p.author = "Ricardo H. Ramirez-Gonzalez"
- p.summary = "Binder of samtools for ruby, on the top of FFI. "
- p.url = "http://github.com/homonecloco/samtools-ruby"
-end
+#Echoe.new("samtools-ruby", version=0.1) do |p|
+# p.author = "Ricardo H. Ramirez-Gonzalez"
+# p.summary = "Binder of samtools for ruby, on the top of FFI. "
+# p.url = "http://github.com/homonecloco/samtools-ruby"
+ # p.platform = Gem
+#end
desc "Basic test"
View
@@ -102,12 +102,38 @@ def load_reference()
end
+ def average_coverage(chromosome, qstart, len)
+ reference = fetch_reference(chromosome, qstart, 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)
+
+ avg_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
+ first.upto(last-1) { |i|
+ coverages[first-i] = 1 + coverages[first-i] }
+ end
+
+ fetch_with_function(chromosome, qstart, qstart+len, avg_proc)
+ p coverages
+ total = 0
+ 0.upto(len-1) { |i| total= total + coverages[i] }
+ avg_cov = total.to_f / len
+ #LibC.free referemce
+ avg_cov
+ end
+
def fetch_reference(chromosome, qstart,qend)
- load_reference unless @fasta_index.nil? || @fasta_index.null?
+ 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
View
@@ -87,6 +87,30 @@ def test_class_binary_read_no_close
assert(true, "BINARY file openend but not closed")
end
+ def test_read_coverage
+ sam = Bio::DB::Sam.new({:bam=>@testBAMFile, :fasta=>@testReference})
+ sam.open
+ File.open( @test_folder +"/ids2.txt", "r") do |file|
+ puts "file opened"
+ file.each_line{|line|
+ fetching = line.split(' ')[0]
+ puts "fetching: " + fetching
+ sam.load_reference
+ seq = sam.fetch_reference(fetching, 0, 16000)
+ # puts seq
+ # puts seq.length
+ als = sam.fetch(fetching, 0, seq.length)
+ # p als
+ if als.length() > 0 then
+ p fetching
+ p als
+ end
+ }
+
+ end
+ sam.close
+ assert(true, "Finish")
+ end
# def test_read_TAM_as_BAM
# begin
# sam = Bio::DB::Sam.new({:bam=>@testTAMFile})
@@ -228,6 +252,17 @@ def test_load_feature
p fs
assert(true, "Loaded as features")
end
+
+ def test_avg_coverage
+ sam = Bio::DB::Sam.new({:fasta=>@testReference, :bam=>@testBAMFile })
+ sam.open
+ cov = sam.average_coverage("chr_1", 60, 30)
+ p "Coverage: " + cov.to_s
+ sam.close
+ assert(true, "Average coverage ran")
+ assert(3 == cov, "The coverage is 3")
+ end
+
end
@@ -0,0 +1 @@
+chr_1

0 comments on commit fd0cd99

Please sign in to comment.