text: #include "global.h"
  #include "aligned_call.h"
  #include "dna.h"
  #include "carrays/carrays.h"
  
  //
  // AlignedCall
  //
  
  void acall_resize(AlignedCall *call, size_t n_lines, size_t n_samples)
  {
    gca_resize(call->lines, call->s_lines, n_lines);
    gca_resize(call->gts, call->s_gts, n_lines*n_samples);
    call->n_lines = n_lines;
    call->n_samples = n_samples;
    call->n_gts = n_lines*n_samples;
  }
  
  void acall_destroy(AlignedCall *call)
  {
    size_t i;
    for(i = 0; i < call->n_lines; i++) strbuf_dealloc(&call->lines[i]);
    free(call->lines);
    free(call->gts);
    strbuf_dealloc(&call->info);
    ctx_free(call);
  }
  
  
  //
  // Decomposer
  //
  
  struct CallDecompStruct
  {
    nw_aligner_t *nw_aligner;
    scoring_t *scoring;
    alignment_t *aln;
    htsFile *vcffh;
    bcf_hdr_t *vcfhdr;
    bcf1_t *v;
    StrBuf sbuf;
    DecomposeStats stats;
  };
  
  CallDecomp* call_decomp_init(htsFile *vcffh, bcf_hdr_t *vcfhdr)
  {
    CallDecomp *dc = ctx_calloc(1, sizeof(CallDecomp));
    dc->nw_aligner = needleman_wunsch_new();
    dc->aln = alignment_create(1024);
    dc->scoring = ctx_calloc(1, sizeof(dc->scoring[0]));
    scoring_system_default(dc->scoring);
    dc->vcffh = vcffh;
    dc->vcfhdr = vcfhdr;
    dc->v = bcf_init();
    strbuf_alloc(&dc->sbuf, 256);
    return dc;
  }
  
  void call_decomp_destroy(CallDecomp *dc)
  {
    alignment_free(dc->aln);
    needleman_wunsch_free(dc->nw_aligner);
    ctx_free(dc->scoring);
    bcf_destroy(dc->v);
    strbuf_dealloc(&dc->sbuf);
    ctx_free(dc);
  }
  
  scoring_t* call_decomp_get_scoring(CallDecomp *dc)
  {
    return dc->scoring;
  }
  
  void call_decomp_cpy_stats(DecomposeStats *stats, const CallDecomp *dc)
  {
    memcpy(stats, &dc->stats, sizeof(*stats));
  }
  
  //
  // Decompose AlignedCall
  //
  
  // Strip indels ('-') from allele and add to string buffer
  static inline void print_vcf_allele(const char *allele, size_t len,
                                      int8_t prev_base, int8_t next_base,
                                      StrBuf *sbuf)
  {
    size_t i;
    if(prev_base > 0) strbuf_append_char(sbuf, char_to_vcf_char(prev_base));
    strbuf_ensure_capacity(sbuf, sbuf->end+len);
    for(i = 0; i < len; i++) {
      if(allele[i] != '-')
        sbuf->b[sbuf->end++] = char_to_vcf_char(allele[i]);
    }
    sbuf->b[sbuf->end] = 0;
    if(next_base > 0) strbuf_append_char(sbuf, char_to_vcf_char(next_base));
  }
  
  // @param vcf_pos is 0-based
  // @param prev_base is -1 if SNP otherwise previous base
  // @param next_base is -1 unless indel at position 0
  static void print_vcf_entry(size_t vcf_pos, int8_t prev_base, int8_t next_base,
                              const char *ref, const char *alt, size_t len,
                              const uint8_t *gts, size_t nsamples,
                              CallDecomp *dc, const AlignedCall *call,
                              size_t max_allele_len)
  {
    dc->stats.nvars++;
  
    StrBuf *sbuf = &dc->sbuf;
    strbuf_reset(sbuf);
  
    // Check actual allele length
    size_t i, alt_bases = 0;
    for(i = 0; i < len; i++) alt_bases += (alt[i] != '-');
    if(alt_bases > max_allele_len) { dc->stats.nallele_too_long++; return; }
  
    // CHROM POS ID REF ALT QUAL FILTER INFO
    strbuf_append_str(sbuf, call->chrom->name.b);
    strbuf_append_char(sbuf, '\t');
    strbuf_append_ulong(sbuf, vcf_pos+1);
    strbuf_append_str(sbuf, "\t.\t");
    print_vcf_allele(ref, len, prev_base, next_base, sbuf);
    strbuf_append_char(sbuf, '\t');
    print_vcf_allele(alt, len, prev_base, next_base, sbuf);
    strbuf_append_str(sbuf, "\t.\tPASS\t");
    strbuf_append_str(sbuf, call->info.b ? call->info.b : ".");
    strbuf_append_str(sbuf, "\tGT");
  
    // Print genotypes
    for(i = 0; i < nsamples; i++) {
      strbuf_append_char(sbuf, '\t');
      strbuf_append_char(sbuf, gts[i] ? '1' : '.');
    }
  
    // fprintf(stderr, " prev_base:%i next_base:%i info:%s\n", prev_base, next_base, call->info.b);
    // fprintf(stderr, "%s [%zu vs %zu]\n", sbuf->b, sbuf->end, strlen(sbuf->b));
  
    kstring_t ks = {.l = sbuf->end, .m = sbuf->size, .s = sbuf->b};
    if(vcf_parse(&ks, dc->vcfhdr, dc->v) != 0)
      die("Cannot construct VCF entry: %s", sbuf->b);
    if(bcf_write(dc->vcffh, dc->vcfhdr, dc->v) != 0)
      die("Cannot write VCF entry [nsamples: %zu vs %zu]", nsamples, (size_t)bcf_hdr_nsamples(dc->vcfhdr));
    // Move back into our string buffer
    sbuf->b = ks.s;
    sbuf->size = ks.m;
  
    dc->stats.nvars_printed++;
  }
  
  // `ref` and `alt` are aligned alleles - should both be same length strings
  // of 'ACGT-'
  // return first mismatch position or -1
  static int align_get_start(const char *ref, const char *alt)
  {
    const char *start = ref;
    while(*ref) {
      if(*ref != *alt) return (ref - start);
      ref++; alt++;
    }
    return -1;
  }
  
  // `ref` and `alt` are aligned alleles - should both be same length strings
  // of 'ACGT-'
  // return first matching position
  static int align_get_end(const char *ref, const char *alt)
  {
    int i = 0;
    while(ref[i] && ref[i] != alt[i]) i++;
    return i;
  }
  
  static size_t align_get_nbases(const char *allele, size_t len)
  {
    size_t i, nbases = 0;
    for(i = 0; i < len; i++)
      if(allele[i] != '-')
        nbases++;
    return nbases;
  }
  
  /**
   * @param ref_pos is 0-based here
   * @param info is extra text to print in the info field of each variant (may be NULL)
   * @param genotypes is strings to print in genotypes columns, of length num_samples.
   *                  It may be NULL.
   * @return number of variants printed
   */
  static void align_biallelic(const char *ref, const char *alt,
                              const read_t *chrom,
                              const uint8_t *gts, size_t nsamples,
                              CallDecomp *dc, const AlignedCall *call,
                              size_t max_allele_len)
  {
    int32_t start, len;
    size_t ref_nbases, alt_nbases, ref_pos = call->start, ref_end, vcf_pos;
    int8_t prev_base, next_base;
    bool is_snp;
  
    // printf("--\n ref: %s\n alt: %s\n", ref, alt);
  
    while((start = align_get_start(ref, alt)) > -1)
    {
      ref_pos += start; // assume ref[i]==alt[i] means ref[i]!='-'
      ref += start;
      alt += start;
      len = align_get_end(ref, alt);
  
      // printf("ref: %.*s\nalt: %.*s\nref_pos: %zu start: %i len %i\n",
      //        len, ref, len, alt, ref_pos, start, len);
  
      ref_nbases = align_get_nbases(ref, len);
      alt_nbases = align_get_nbases(alt, len);
      is_snp = (ref_nbases == 1 && alt_nbases == 1);
      ref_end = ref_pos+ref_nbases;
      vcf_pos = ref_pos; // copy in case we need left padding base
  
      // If one allele is going to be empty, we need a padding base
      // If ref_pos == 0, add extra base to end instead
      prev_base = next_base = -1;
      if(!is_snp) {
        if(ref_pos > 0) prev_base = chrom->seq.b[--vcf_pos];
        else if(ref_end < chrom->seq.end) next_base = chrom->seq.b[ref_end];
      }
  
      if(is_snp || prev_base > 0 || next_base > 0) {
        print_vcf_entry(vcf_pos, prev_base, next_base, ref, alt, len,
                        gts, nsamples, dc, call, max_allele_len);
      }
  
      ref_pos += ref_nbases;
      ref += len;
      alt += len;
    }
  }
  
  void acall_decompose(CallDecomp *dc, const AlignedCall *call,
                       size_t max_line_len, size_t max_allele_len)
  {
    dc->stats.ncalls++;
    if(call->chrom == NULL) { return; }
    dc->stats.ncalls_mapped++;
  
    const read_t *chrom = call->chrom;
    const char *ref_allele = chrom->seq.b + call->start;
    size_t i, ref_len = call->end - call->start;
    const StrBuf *alt;
  
    ctx_assert2(call->start <= call->end, "%u .. %u", call->start, call->end);
  
    if(ref_len > max_line_len) {
      dc->stats.ncalls_ref_allele_too_long++;
      return; // can't align
    }
  
    dc->stats.nlines += call->n_lines;
  
    // printf("chr:%s %u - %u\n", call->chrom->name.b, call->start, call->end);
  
    for(i = 0; i < call->n_lines; i++)
    {
      alt = &call->lines[i];
      ctx_assert(strlen(alt->b) == alt->end);
  
      // Quick check if sequence too long or are matching
      if(alt->end > max_line_len) {
        dc->stats.nlines_too_long++;
      } else if(ref_len == alt->end && strncasecmp(ref_allele, alt->b, ref_len) == 0) {
        dc->stats.nlines_match_ref++;
      } else {
        // printf("REF: '%*.s' [%zu]\n", (int)ref_len, ref_allele, ref_len);
        // printf("ALT: '%*.s' [%zu]\n", (int)alt->end, alt->b, alt->end);
  
        needleman_wunsch_align2(ref_allele, alt->b, ref_len, alt->end,
                                dc->scoring, dc->nw_aligner, dc->aln);
  
        // printf("ALNA: %s\n", dc->aln->result_a);
        // printf("ALNB: %s\n", dc->aln->result_b);
  
        align_biallelic(dc->aln->result_a, dc->aln->result_b, chrom,
                        call->gts+i*call->n_samples, call->n_samples,
                        dc, call, max_allele_len);
        dc->stats.nlines_mapped++;
      }
    }
  }
  