#!/bin/bash -Ceuo pipefail
snv=$(bcftools view -H -f PASS -v snps   germline_control_2.haplotyper_rf.gt_corrected.vcf.gz | wc -l)
indel=$(bcftools view -H -f PASS -v indels germline_control_2.haplotyper_rf.gt_corrected.vcf.gz | wc -l)

{
    printf "snv_count\t%s\n"   "$snv"
    printf "indel_count\t%s\n" "$indel"
} > germline_control_2.snv_indel.tsv

bcftools norm -m - -f Homo_sapiens_assembly38.fasta -Oz -o giab.norm.vcf.gz na12878_hardcore.vcf.gz
bcftools index -t giab.norm.vcf.gz

bcftools view -f PASS -Oz -o calls.filt.vcf.gz germline_control_2.haplotyper_rf.gt_corrected.vcf.gz
bcftools norm -m - -f Homo_sapiens_assembly38.fasta -Oz -o calls.norm.vcf.gz calls.filt.vcf.gz
bcftools index -t calls.norm.vcf.gz

bcftools isec -p isec_out -Oz calls.norm.vcf.gz giab.norm.vcf.gz

giab_total=$(bcftools view -H giab.norm.vcf.gz      | wc -l)
giab_recov=$(bcftools view -H isec_out/0003.vcf.gz  | wc -l)
giab_missed=$(bcftools query -f '%ID|%CHROM:%POS:%REF:%ALT\n' isec_out/0001.vcf.gz \
    | awk -F'|' '{ if ($1 != "" && $1 != ".") print $1; else print $2 }' \
    | paste -sd';' -)
: ${giab_missed:=NONE}

{
    printf "giab_truth_total\t%s\n"     "$giab_total"
    printf "giab_recovered_count\t%s\n" "$giab_recov"
    printf "giab_missed_variants\t%s\n" "$giab_missed"
} >> germline_control_2.snv_indel.tsv


cat <<-END_VERSIONS > versions.yml
"DAQ:CONTROL_METRICS:GERMLINE_METRICS:GERMLINE_SNV_INDEL":
    bcftools: $(bcftools --version | head -n1 | awk '{print $2}')
END_VERSIONS