Add auxiliary tag manipulation examples and tests

kojix2 · kojix2 · commit 1f50b350d81e · 2025-11-27T10:17:58.000Z
diff --git a/TUTORIAL.md b/TUTORIAL.md
@@ -254,6 +254,51 @@ in.close
 out.close
 ```
 
+Writing and modifying auxiliary tags
+
+```ruby
+# Reading auxiliary tags
+bam = HTS::Bam.open("input.bam")
+record = bam.first
+aux = record.aux
+
+# Read tags
+alignment_score = aux["AS"]           # Auto-detect type
+mc_cigar = aux.get_string("MC")       # Type-specific getter
+edit_distance = aux.get_int("NM")     # Type-specific getter
+
+# Writing/updating auxiliary tags
+in_bam = HTS::Bam.open("input.bam")
+out_bam = HTS::Bam.open("output.bam", "wb")
+out_bam.write_header(in_bam.header)
+
+in_bam.each do |record|
+  aux = record.aux
+  
+  # Update or add tags using type-specific methods
+  aux.update_int("AS", 100)                    # Integer tag
+  aux.update_float("ZQ", 0.95)                 # Float tag
+  aux.update_string("RG", "sample1")           # String tag
+  aux.update_array("BC", [25, 30, 28, 32])     # Array tag
+  
+  # Or use the []= operator (auto-detects type)
+  aux["NM"] = 2                                # Integer
+  aux["ZS"] = "modified"                       # String
+  aux["ZF"] = 3.14                             # Float
+  aux["ZA"] = [1, 2, 3, 4]                     # Array
+  
+  # Check if tag exists
+  if aux.key?("XS")
+    aux.delete("XS")  # Delete tag
+  end
+  
+  out_bam.write(record)
+end
+
+in_bam.close
+out_bam.close
+```
+
 Create index
 
 ```ruby
diff --git a/examples/bam-write.rb b/examples/bam-write.rb
@@ -0,0 +1,119 @@
+#!/usr/bin/env ruby
+# frozen_string_literal: true
+
+# Example: Writing BAM files
+#
+# This example demonstrates how to:
+# - Read BAM files and write modified records to a new BAM file
+# - Modify record fields (MAPQ, flags, positions, etc.)
+# - Add, update, and delete auxiliary tags
+# - Filter records based on criteria
+#
+# Usage: bam-write.rb <input.bam> <output.bam>
+
+require "htslib"
+
+def main
+  if ARGV.size < 2
+    warn "Usage: #{$PROGRAM_NAME} <input.bam> <output.bam>"
+    warn ""
+    warn "This example reads a BAM file, modifies records, and writes to a new file."
+    warn ""
+    warn "Examples:"
+    warn "  #{$PROGRAM_NAME} input.bam output.bam"
+    exit 1
+  end
+
+  input_path = ARGV[0]
+  output_path = ARGV[1]
+
+  # Open input BAM file
+  input_bam = HTS::Bam.new(input_path)
+  header = input_bam.header
+
+  # Open output BAM file for writing
+  output_bam = HTS::Bam.new(output_path, "wb")
+  output_bam.write_header(header)
+
+  written_count = 0
+  filtered_count = 0
+
+  # Process each record
+  input_bam.each do |record|
+    # Example 1: Filter records by mapping quality
+    # Skip records with low mapping quality
+    if record.mapq < 10
+      filtered_count += 1
+      next
+    end
+
+    # Example 2: Modify record fields
+    # Adjust mapping quality for high-confidence alignments
+    if record.mapq > 50
+      record.mapq = 60 # Cap at 60
+    end
+
+    # Example 3: Modify flags
+    # Mark duplicates (example logic)
+    if record.pos % 100 == 0 # Arbitrary example
+      record.flag
+      # This is just an example - real duplicate marking is more complex
+      # flag |= 1024  # Set duplicate flag
+    end
+
+    # Example 4: Work with auxiliary tags
+    aux = record.aux
+
+    # Add or update tags using type-specific methods
+    aux.update_int("AS", record.mapq * 2)  # Alignment score based on MAPQ
+    aux.update_string("PG", "bam-write")   # Program name
+
+    # Or use the []= operator (auto-detects type)
+    aux["NM"] = 0              # Reset edit distance (example)
+    aux["ZP"] = "processed"    # Custom processing flag
+
+    # Add custom quality metrics
+    if record.mapq > 30
+      aux.update_float("ZQ", 0.95)  # High quality score
+      aux["HQ"] = 1                 # High quality flag
+    else
+      aux.update_float("ZQ", 0.50)  # Medium quality score
+    end
+
+    # Delete tags if needed (example: remove XS tag)
+    aux.delete("XS") if aux.key?("XS")
+
+    # Example 5: Add array tags
+    # Add custom base quality distribution (mock data)
+    if record.seq && record.seq.length > 0
+      # Count ACGT bases (simplified example)
+      seq = record.seq
+      base_counts = [
+        seq.count("A"),
+        seq.count("C"),
+        seq.count("G"),
+        seq.count("T")
+      ]
+      aux.update_array("BC", base_counts)
+    end
+
+    # Write modified record to output
+    output_bam.write(record)
+    written_count += 1
+
+    # Print progress every 10000 records
+    warn "Processed: #{written_count} written, #{filtered_count} filtered..." if (written_count % 10_000).zero?
+  end
+
+  # Clean up
+  input_bam.close
+  output_bam.close
+
+  puts "Done!"
+  puts "  Records written: #{written_count}"
+  puts "  Records filtered: #{filtered_count}"
+  puts "  Total processed: #{written_count + filtered_count}"
+  puts "  Output: #{output_path}"
+end
+
+main if $PROGRAM_NAME == __FILE__
diff --git a/lib/hts/bam/auxi.rb b/lib/hts/bam/auxi.rb
@@ -55,6 +55,124 @@ def [](key)
         get(key)
       end
 
+      # Set auxiliary tag value (auto-detects type from value)
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @param value [Integer, Float, String, Array] tag value
+      def []=(key, value)
+        case value
+        when Integer
+          update_int(key, value)
+        when Float
+          update_float(key, value)
+        when String
+          update_string(key, value)
+        when Array
+          update_array(key, value)
+        else
+          raise ArgumentError, "Unsupported type: #{value.class}"
+        end
+      end
+
+      # Update or add an integer tag
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @param value [Integer] integer value
+      def update_int(key, value)
+        ret = LibHTS.bam_aux_update_int(@record.struct, key, value.to_i)
+        raise "Failed to update integer tag '#{key}': errno #{FFI.errno}" if ret < 0
+
+        value
+      end
+
+      # Update or add a floating-point tag
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @param value [Float] floating-point value
+      def update_float(key, value)
+        ret = LibHTS.bam_aux_update_float(@record.struct, key, value.to_f)
+        raise "Failed to update float tag '#{key}': errno #{FFI.errno}" if ret < 0
+
+        value
+      end
+
+      # Update or add a string tag
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @param value [String] string value
+      def update_string(key, value)
+        ret = LibHTS.bam_aux_update_str(@record.struct, key, -1, value.to_s)
+        raise "Failed to update string tag '#{key}': errno #{FFI.errno}" if ret < 0
+
+        value
+      end
+
+      # Update or add an array tag
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @param value [Array] array of integers or floats
+      # @param type [String, nil] element type ('c', 'C', 's', 'S', 'i', 'I', 'f'). Auto-detected if nil.
+      def update_array(key, value, type: nil)
+        raise ArgumentError, "Array cannot be empty" if value.empty?
+
+        # Auto-detect type if not specified
+        if type.nil?
+          if value.all? { |v| v.is_a?(Integer) }
+            # Use 'i' for signed 32-bit integers by default
+            type = "i"
+          elsif value.all? { |v| v.is_a?(Float) || v.is_a?(Integer) }
+            type = "f"
+          else
+            raise ArgumentError, "Array must contain only integers or floats"
+          end
+        end
+
+        # Convert array to appropriate C type
+        case type
+        when "c", "C", "s", "S", "i", "I"
+          # Integer types
+          ptr = FFI::MemoryPointer.new(:int32, value.size)
+          ptr.write_array_of_int32(value.map(&:to_i))
+          ret = LibHTS.bam_aux_update_array(@record.struct, key, type.ord, value.size, ptr)
+        when "f"
+          # Float type
+          ptr = FFI::MemoryPointer.new(:float, value.size)
+          ptr.write_array_of_float(value.map(&:to_f))
+          ret = LibHTS.bam_aux_update_array(@record.struct, key, type.ord, value.size, ptr)
+        else
+          raise ArgumentError, "Invalid array type: #{type}"
+        end
+
+        raise "Failed to update array tag '#{key}': errno #{FFI.errno}" if ret < 0
+
+        value
+      end
+
+      # Delete an auxiliary tag
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @return [Boolean] true if tag was deleted, false if tag was not found
+      def delete(key)
+        aux_ptr = LibHTS.bam_aux_get(@record.struct, key)
+        return false if aux_ptr.null?
+
+        ret = LibHTS.bam_aux_del(@record.struct, aux_ptr)
+        raise "Failed to delete tag '#{key}': errno #{FFI.errno}" if ret < 0
+
+        true
+      end
+
+      # Check if a tag exists
+      # For compatibility with HTS.cr.
+      # @param key [String] tag name (2 characters)
+      # @return [Boolean] true if tag exists
+      def key?(key)
+        aux_ptr = LibHTS.bam_aux_get(@record.struct, key)
+        !aux_ptr.null?
+      end
+
+      alias include? key?
+
       def first
         aux_ptr = first_pointer
         return nil if aux_ptr.null?
diff --git a/test/bam_test.rb b/test/bam_test.rb
