import hts

# open a bam and look for the index.
var b:Bam
open(b, "tests/HG02002.bam", index=true)

for record in b:
  if record.qual > 10u:
    echo record.chrom, record.start, record.stop

# regional queries:
for record in b.query("6", 30816675, 32816675):
  if record.flag.proper_pair and record.flag.reverse:
    # cigar is an iterable of operations:
    for op in record.cigar:
      # $op gives the string repr of the operation, e.g. '151M'
      echo $op, op.consumes.reference, op.consumes.query

    # tags are pulled with `aux`
    var mismatches = rec.aux("NM")
    if mismatches != nil and mismatches.asInt.get < 3:
      var rg = rec.aux("RG")
      echo rg.asString

# cram requires an fasta to decode:
var cram:Bam
open(cram, "/tmp/t.cram", fai="/data/human/g1k_v37_decoy.fa")
for record in cram:
  # now record is same as from bam above
  echo record.qname, record.isize

import hts

var tsamples = @["101976-101976", "100920-100920", "100231-100231", "100232-100232", "100919-100919"]

# VCF and BCF supported
var v:VCF
assert open(v, "tests/test.bcf", samples=tsamples)

var afs = new_seq[float32](5) # size doesn't matter. this will be re-sized as needed
var acs = new_seq[int32](5) # size doesn't matter. this will be re-sized as needed
var csq = new_string_of_cap(20)
for rec in v:
  echo rec, " qual:", rec.QUAL, " filter:", rec.FILTER
  var info = rec.info
  # accessing stuff from the INFO field is meant to be as fast as possible, allowing
  # the user to re-use memory as needed.
  info.strings("CSQ", csq)
  info.ints("AC", acs)
  info.floats("AF", afs)
  echo acs, afs, csq, info.has_flag("IN_EXAC")

  # accessing format fields is similar
  var format = ref.format
  var dps = new_seq[int32](len(v.samples))
  assert format.ints("DP", dps)
  echo dps

echo v.samples

# open a VCF for writing
var wtr:VCF
open(wtr, "tests/outv.vcf", mode="w")
# set and write the header.
wtr.header = v.header
assert wtr.write_header()

# regional queries look for index. works for VCF and BCF
for rec in v.query("1:15600-18250"):
  echo rec.CHROM, ":", $rec.POS
  # adjust some values in the INFO
  var val = 22.3
  check rec.info.set("VQSLOD", val) == Status.OK
  assert wtr.write_variant(rec)

  check rec.info.delete("CSQ") == Status.OK

wtr.close()

import hts
# arguments are 1-based seq-col, start-col, end-col and whether the intervals are 0-based.
var bx = wopen_bgzi("ti.txt.gz", 1, 2, 3, true)
# some duplication of args to avoid re-parsing. args are line, chrom, start, end
check bx.write_interval("a\t4\t10", "a", 4, 10) > 0
check bx.write_interval("b\t2\t20", "b", 2, 20) > 0
check bx.close() == 0