Enhance BCF format handling with new decoding methods for integers and floats, and add tests for high-level value retrieval and raw buffer preservation

Author: kojix2

lib/hts/bcf/format.rb

@@ -12,6 +12,23 @@ def initialize(record)
       # which provides methods like `get_int`, `get_float`, etc.
       # I think they are better than `fetch_int`` and `fetch_float`.
       def get(key, type = nil)
+        return get_raw(key, type) unless type.nil?
+
+        return decode_genotypes if key == "GT"
+
+        case header_format_type(key)
+        when :int
+          decode_integer_values(key)
+        when :float
+          decode_float_values(key)
+        when :flag
+          raise_unsupported_format_flag(key)
+        when :string
+          get_string_values(key)
+        end
+      end
+
+      def get_raw(key, type = nil)
         # The GT FORMAT field is special in that it is marked as a string in the header,
         # but it is actually encoded as an integer.
         type = if type.nil?
@@ -39,22 +56,22 @@ def get(key, type = nil)
 
       # For compatibility with HTS.cr.
       def get_int(key)
-        get(key, :int)
+        get_raw(key, :int)
       end
 
       # For compatibility with HTS.cr.
       def get_float(key)
-        get(key, :float)
+        get_raw(key, :float)
       end
 
       # For compatibility with HTS.cr.
       def get_flag(key)
-        get(key, :flag)
+        get_raw(key, :flag)
       end
 
       # For compatibility with HTS.cr.
       def get_string(key)
-        get(key, :string)
+        get_raw(key, :string)
       end
 
       # For compatibility with HTS.cr.
@@ -201,6 +218,38 @@ def get_string_values(key)
         end
       end
 
+      def decode_integer_values(key)
+        values = get_raw(key, :int)
+        return nil unless values
+
+        sample_values = split_sample_values(values)
+        if scalar_format?(key)
+          sample_values.map do |values_per_sample|
+            map_integer_missing_value(trim_integer_vector_end(values_per_sample).first)
+          end
+        else
+          sample_values.map do |values_per_sample|
+            map_integer_missing(trim_integer_vector_end(values_per_sample))
+          end
+        end
+      end
+
+      def decode_float_values(key)
+        values = get_float_words(key)
+        return nil unless values
+
+        sample_values = split_sample_values(values)
+        if scalar_format?(key)
+          sample_values.map do |values_per_sample|
+            decode_float_word(trim_float_vector_end(values_per_sample).first)
+          end
+        else
+          sample_values.map do |values_per_sample|
+            map_float_words(trim_float_vector_end(values_per_sample))
+          end
+        end
+      end
+
       def decode_genotypes
         genotypes = get_genotypes
         return nil unless genotypes
@@ -245,6 +294,38 @@ def trim_genotype_vector_end(values)
         values[0, end_index]
       end
 
+      def trim_integer_vector_end(values)
+        end_index = values.index { |value| value == LibHTS.bcf_int32_vector_end } || values.size
+        values[0, end_index]
+      end
+
+      def trim_float_vector_end(values)
+        end_index = values.index(0x7f80_0002) || values.size
+        values[0, end_index]
+      end
+
+      def map_integer_missing(values)
+        values.map { |value| map_integer_missing_value(value) }
+      end
+
+      def map_integer_missing_value(value)
+        value == LibHTS.bcf_int32_missing ? nil : value
+      end
+
+      def map_float_words(values)
+        values.map { |value| decode_float_word(value) }
+      end
+
+      def decode_float_word(value)
+        return nil if value == 0x7f80_0001
+
+        [value].pack("V").unpack1("e")
+      end
+
+      def get_float_words(key)
+        get_numeric_values(key, LibHTS::BCF_HT_REAL, "float") { |dst, len| dst.get_array_of_uint32(0, len) }
+      end
+
       def normalize_format_rc(rc, key, expected_type)
         case rc
         when -1, -3
@@ -364,6 +445,17 @@ def get_fmt_type(qname)
         nil
       end
 
+      def scalar_format?(key)
+        header_format_number(key) == 1
+      end
+
+      def header_format_number(key)
+        id = LibHTS.bcf_hdr_id2int(@record.header.struct, LibHTS::BCF_DT_ID, key)
+        return nil unless LibHTS.bcf_hdr_idinfo_exists(@record.header.struct, LibHTS::BCF_HL_FMT, id)
+
+        LibHTS.bcf_hdr_id2number(@record.header.struct, LibHTS::BCF_HL_FMT, id)
+      end
+
       def header_format_type_code(key)
         id = LibHTS.bcf_hdr_id2int(@record.header.struct, LibHTS::BCF_DT_ID, key)
         return nil unless LibHTS.bcf_hdr_idinfo_exists(@record.header.struct, LibHTS::BCF_HL_FMT, id)
@@ -388,6 +480,7 @@ def ht_type_to_sym(t)
       def int32_range?(value)
         value >= -2_147_483_648 && value <= 2_147_483_647
       end
+
     end
   end
 end

test/bcf/format_test.rb

@@ -36,6 +36,69 @@ def with_temp_character_format_vcf
     end
   end
 
+  def with_temp_bcf
+    Tempfile.create(["format_test", ".bcf"]) do |file|
+      path = file.path
+      file.close
+
+      header = HTS::Bcf::Header.new
+      header.set_version("VCFv4.3")
+      header.append("##contig=<ID=1,length=100>")
+      header.append('##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">')
+      header.append('##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Phred likelihoods">')
+      header.append('##FORMAT=<ID=MISSI,Number=1,Type=Integer,Description="defined but absent integer">')
+      header.append('##FORMAT=<ID=MISSF,Number=1,Type=Float,Description="defined but absent float">')
+      header.append('##FORMAT=<ID=IV,Number=.,Type=Integer,Description="integer with sentinels">')
+      header.append('##FORMAT=<ID=FV,Number=.,Type=Float,Description="float with sentinels">')
+      header.add_sample("S1", sync: false)
+      header.add_sample("S2", sync: true)
+
+      HTS::Bcf.open(path, "wb") do |bcf|
+        bcf.write_header(header)
+
+        record = HTS::Bcf::Record.new(header)
+        record.rid = HTS::LibHTS.bcf_hdr_name2id(header.struct, "1")
+        record.pos = 9
+
+        rc = HTS::LibHTS.bcf_update_alleles_str(header.struct, record.struct, "A,C")
+        raise "bcf_update_alleles_str failed (rc=#{rc})" if rc < 0
+
+        genotypes = [
+          HTS::LibHTS.bcf_gt_unphased(0),
+          HTS::LibHTS.bcf_gt_unphased(1),
+          HTS::LibHTS.bcf_gt_unphased(1),
+          HTS::LibHTS.bcf_gt_unphased(1)
+        ]
+        genotype_ptr = FFI::MemoryPointer.new(:int32, genotypes.size)
+        genotype_ptr.write_array_of_int32(genotypes)
+        rc = HTS::LibHTS.bcf_update_genotypes(header.struct, record.struct, genotype_ptr, genotypes.size)
+        raise "bcf_update_genotypes failed (rc=#{rc})" if rc < 0
+
+        likelihoods = [10, 20, 30, 40, 50, 60]
+        likelihood_ptr = FFI::MemoryPointer.new(:int32, likelihoods.size)
+        likelihood_ptr.write_array_of_int32(likelihoods)
+        rc = HTS::LibHTS.bcf_update_format_int32(header.struct, record.struct, "PL", likelihood_ptr, likelihoods.size)
+        raise "bcf_update_format_int32 failed (rc=#{rc})" if rc < 0
+
+        int_with_sentinels = [10, HTS::LibHTS.bcf_int32_vector_end, HTS::LibHTS.bcf_int32_missing, HTS::LibHTS.bcf_int32_vector_end]
+        int_ptr = FFI::MemoryPointer.new(:int32, int_with_sentinels.size)
+        int_ptr.write_array_of_int32(int_with_sentinels)
+        rc = HTS::LibHTS.bcf_update_format_int32(header.struct, record.struct, "IV", int_ptr, int_with_sentinels.size)
+        raise "bcf_update_format_int32 failed for IV (rc=#{rc})" if rc < 0
+
+        float_words = [0x3fc0_0000, 0x7f80_0002, 0x7f80_0001, 0x7f80_0002]
+        float_ptr = FFI::MemoryPointer.new(:uint32, float_words.size)
+        float_ptr.write_array_of_uint32(float_words)
+        rc = HTS::LibHTS.bcf_update_format_float(header.struct, record.struct, "FV", float_ptr, float_words.size)
+        raise "bcf_update_format_float failed for FV (rc=#{rc})" if rc < 0
+
+        bcf.write(record)
+      end
+
+      yield path
+    end
+  end
+
   def with_temp_gt_vcf
     Tempfile.create(["format_gt", ".vcf"]) do |file|
       file.write <<~VCF
@@ -112,13 +175,15 @@ def test_get_like_crystal
   def test_get_without_type
     assert_equal [409, 409], @fmt.get("GQ")
     assert_equal [35, 35], @fmt.get("DP")
-    assert_equal [-20.0, -5.0, -20.0, -20.0, -5.0, -20.0], @fmt.get("GL")
+    assert_equal [[-20.0, -5.0, -20.0], [-20.0, -5.0, -20.0]], @fmt.get("GL")
+    assert_equal ["0/1", "0/1"], @fmt.get("GT")
   end
 
   def test_get_square_brackets
     assert_equal [409, 409], @fmt["GQ"]
     assert_equal [35, 35], @fmt["DP"]
-    assert_equal [-20.0, -5.0, -20.0, -20.0, -5.0, -20.0], @fmt["GL"]
+    assert_equal [[-20.0, -5.0, -20.0], [-20.0, -5.0, -20.0]], @fmt["GL"]
+    assert_equal ["0/1", "0/1"], @fmt["GT"]
   end
 
   def test_get_unknown_key
@@ -189,13 +254,55 @@ def test_fields
 
   def test_to_h
     assert_equal(
-      # FIXME: Maybe GT should be string?
-      { "GT" => [2, 4, 2, 4], "GQ" => [409, 409], "DP" => [35, 35],
-        "GL" => [-20.0, -5.0, -20.0, -20.0, -5.0, -20.0] },
+      { "GT" => ["0/1", "0/1"], "GQ" => [409, 409], "DP" => [35, 35],
+        "GL" => [[-20.0, -5.0, -20.0], [-20.0, -5.0, -20.0]] },
       @fmt.to_h
     )
   end
 
+  def test_get_high_level_shapes_values_by_sample
+    with_temp_bcf do |path|
+      HTS::Bcf.open(path) do |bcf|
+        format = bcf.first.format
+
+        assert_equal ["0/1", "1/1"], format.get("GT")
+        assert_equal [[10, 20, 30], [40, 50, 60]], format.get("PL")
+        assert_equal [[10], [nil]], format.get("IV")
+
+        floats = format.get("FV")
+        assert_equal 2, floats.size
+        assert_equal [1.5], floats[0]
+        assert_equal [nil], floats[1]
+        assert_nil format.get("MISSI")
+        assert_nil format.get("MISSF")
+      end
+    end
+  end
+
+  def test_get_raw_preserves_flat_numeric_buffers
+    with_temp_bcf do |path|
+      HTS::Bcf.open(path) do |bcf|
+        format = bcf.first.format
+
+        assert_equal [10, 20, 30, 40, 50, 60], format.get_raw("PL")
+
+        ints = format.get_raw("IV")
+        assert_equal 4, ints.size
+        assert_equal 10, ints[0]
+        assert_equal HTS::LibHTS.bcf_int32_vector_end, ints[1]
+        assert_equal HTS::LibHTS.bcf_int32_missing, ints[2]
+        assert_equal HTS::LibHTS.bcf_int32_vector_end, ints[3]
+
+        floats = format.get_raw("FV")
+        assert_equal 4, floats.size
+        assert_in_delta 1.5, floats[0], 0.001
+        assert_predicate floats[1], :nan?
+        assert_predicate floats[2], :nan?
+        assert_predicate floats[3], :nan?
+      end
+    end
+  end
+
   def test_update_methods_round_trip
     with_temp_format_source_vcf do |source_path|
       with_temp_output_path do |output_path|