#!/bin/bash

##Run from inside needLR_local

##Load environment with 
  #truvari v4.2.2 
  #bedtools v2.31.1 
  #bcftools v1.19 

##Default Truvari parameters
truvari_k="common"
truvari_s=50
truvari_S=10000000
truvari_B=50
truvari_r=2000
truvari_p=0
truvari_P=0.2
truvari_O=0.2

while getopts g:t:l:u:n:s:b:r:m:p: flag
do
    case "${flag}" in
        g) REF_GENOME=${OPTARG};;
        #t) TRUVARI_MERGED_CTRL_SET=${OPTARG};;
        #l) QUERY_VCFS=${OPTARG};;
        u) MULTISAMPLE_VCF=${OPTARG};;
        n) CTRL_SAMPLE_NAMES=${OPTARG};;
        s) TOTAL_SAMPLES=${OPTARG};;
        #b) PROBAND_VCF=${OPTARG};;
        #r) PARENTAL_VCF=${OPTARG};;
        #m) MATERNAL_VCF=${OPTARG};;
        #p) PATERNAL_VCF=${OPTARG};;
    esac
done

##This will be the basename for the output files
QUERY_FILE_NAME=$(basename "$MULTISAMPLE_VCF" .vcf.gz)_needLR_3.4_anno

##This is the output directory specific to each query VCF
UNSOLVED_DIR="needLR_output/$QUERY_FILE_NAME"
mkdir $UNSOLVED_DIR

##These are the bed files used for annotating the SVs
GENES="backend_files/bed_files/PROTEIN_CODING_GENE_gencode.v45.annotation_1kb_slop.bed"
CODING_REGIONS="backend_files/bed_files/ENSEMBL_CANONICAL_EXON_in_PROTEIN_CODING_GENE_gencode.v45.annotation.bed"
OMIM_GENES="backend_files/bed_files/OMIM_gene_phen_hg38.bed"
CENTROMERES="backend_files/bed_files/hg38_centromeres_endtoend.bed"
PERICENTROMERES="backend_files/bed_files/hg38_pericentromeres_5Mb.bed"
TELOMERES="backend_files/bed_files/hg38_telomeres_5Mb.bed"
DEFRABB_HICONF="backend_files/bed_files/GRCh38_HG002-T2TQ100-V1.0_stvar.benchmark.bed"
STR="backend_files/bed_files/STR_original_motifs.set148.bed"
VNTR="backend_files/bed_files/VNTR_original_motifs.set148.bed"
SEGDUPS="backend_files/bed_files/SEGDUPS_GIAB_v3.3.bed"
REPEAT_MASKER="backend_files/bed_files/UCSC_hg38_Repeats_RepeatMasker.bed"
GAPS="backend_files/bed_files/UCSC_hg38_Mapping_and_Sequencing_Gap.bed"

##This parses the Truvari VCF output for relevant information
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%INFO/SVLEN\t%INFO/SVTYPE[\t%GT]\n' "$MULTISAMPLE_VCF" > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_parsed_temp.txt

##This removes BNDs from the analysis (too large to efficiently annotate)
cat "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_parsed_temp.txt | grep -v BND > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp1.txt

##This parses for SVs <=10Mb
awk '($5 > -10000000 && $5 < 10000000) { print }' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp1.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp2.txt

##This creates a SUPP_VEC (like in Jasmine)

#$1 - Chr
#$2 - Pos
#$3 - REF
#$4 - ALT
#$5 - SV len
#$6 - SV type
# #$7-? - CTRl genotypes

TEMPA=$(( TOTAL_SAMPLES + 6 ))

awk -v tempa="$TEMPA" 'BEGIN {OFS="\t"} {
    result = "";
    for (i=7; i<=tempa; i++) {
        if ($i ~ /1/) {
            result = result "1";
        } else {
            result = result "0";
        }
    }
    print $0, result
}' "$UNSOLVED_DIR/$QUERY_FILE_NAME"_temp2.txt \
> "$UNSOLVED_DIR/$QUERY_FILE_NAME"_temp3.txt

##Rearrange columns so SUPP_VEC is $8
##START_COL = first column of control genotypes
START_COL=7
END_COL=$(( TOTAL_SAMPLES + 6 ))
SUPP_VECC=$(( TOTAL_SAMPLES + 7 ))
STEP=1

awk -v start="$START_COL" -v end="$END_COL" -v step="$STEP" -v suppvecc="$SUPP_VECC" 'BEGIN {OFS="\t"}
{printf("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t",
         $1, $2, "temp", $3, $4, $5, $6, $(suppvecc))

  # Print the range of columns in a loop
  for (i = start; i <= end; i += step) {
    printf("%s\t", $i)
  }
  printf("\n")
}
' "$UNSOLVED_DIR/$QUERY_FILE_NAME"_temp3.txt \
> "$UNSOLVED_DIR/$QUERY_FILE_NAME"_temp4.txt


##This calculates the population frequencies, allele frequencies, and allele types (for HWE)

TEMPALL=$((TOTAL_SAMPLES + 8))

awk -v total_samples="$TOTAL_SAMPLES" -v tempall="$TEMPALL" 'BEGIN {OFS="\t"} {
   substring_ALL = substr($8, 1, total_samples)
   
   ## Count the number of 1s in the substring
   Pop_ALL = gsub("1", "", substring_ALL)

   Allele_ALL = 0;
   GT_homWT = 0;
   GT_het = 0;
   GT_homVAR = 0;
  
   for (i = 9; i <= tempall; i++) {
      n = split($i, arr, "");
       for (j = 1; j <=n; j++) {
          if (arr[j] == "1") {
           Allele_ALL++;
       }
     }
   }
   for (i = 9; i <= tempall; i++) {
     if ($i == "./." || $i == "0/0") {
           GT_homWT++;
       }
     }
   for (i = 9; i <= tempall; i++) {
     if ($i == "0/1" || $i == "1/0") {
           GT_het++;
       }
     }
   for (i = 9; i <= tempall; i++) {
     if ($i == "1/1") {
          GT_homVAR++;
       }
     }

   print $1, $2, $3, $4, $5, $6, $7, $8, Pop_ALL, Allele_ALL, GT_homWT, GT_het, GT_homVAR
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp4.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp5.txt

##This matches the SUPP_VEC string to the sample list (and lists which 1KGP samples have the SV)
mapfile -t names < "$CTRL_SAMPLE_NAMES"

awk -v names="$(printf "%s\n" "${names[@]}")" '
BEGIN {
    split(names, nameArray, "\n");
    OFS="\t";
}
{
    # Get the binary string from column 8
    binary_string = $8;
    
    # Check if the binary string is all zeros
    if (binary_string ~ /^0+$/) {
        names_string = "query_only";
    } else {
        # Initialize the names string
        names_string = "";
        
        # Loop over each character in the binary string
        for (i = 1; i <= length(binary_string); i++) {
            if (substr(binary_string, i, 1) == "1") {
                if (names_string == "") {
                    names_string = nameArray[i];
                } else {
                    names_string = names_string "_" nameArray[i];
                }
            }
        }
    }
    
    # Replace column 8 with the names string
    $8 = names_string;
    
    # Print the modified row
    print $0;
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp5.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp6.txt

#$11 - GT_homWT 
#$12 - GT_het 
#$13 - GT_homVAR

# This calculates p for HWE - WT alle freq ((2*homWT) + het)/(2* (homWT + het + homVar))
# Checks for demoninator of zero (if SV is counted but genotype is ./. for a 1KGP sample)

awk '{
  if ($11 + $12 + $13 != 0) {
    result = ((2 * $11) + $12) / (2 * ($11 + $12 + $13));
  } else {
    result = 0;
  }
  print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, result
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp6.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp7.txt

##This calculates q for HWE - (1-p)
awk '{print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, (1-$14)}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp7.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp8.txt

##This calculates np^2, 2npq, and nq^2 for HWE

TOTAL_HAPS=$((TOTAL_SAMPLES * 2))

awk -v total_haps="$TOTAL_HAPS" -v total_samples="$TOTAL_SAMPLES" '{print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, $15, (total_samples*$14*$14), (total_haps*$14*$15), (total_samples*$15*$15)}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp8.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp9.txt

#This calculates chi-squared value for HWE. (Can't divide by zero, so if expected is zero, chi-sq=NULL) 
awk '{
   if ($16 > 0 && $17 > 0 && $18 > 0) {
       result = ((($11 - $16) * ($11 - $16)) / $16) + ((($12 - $17) * ($12 - $17)) / $17) + ((($13 - $18) * ($13 - $18)) / $18);
       print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, $15, $16, $17, $18, result;
   } else {
       print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, $15, $16, $17, $18, "NULL";
   }
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp9.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp10.txt


## $9 - Pop_ALL
## $10 - Allele_ALL
## $11 - GT_homWT
## $12 - GT_het
## $13 - GT_homVAR

TOTAL_HAPS=$((TOTAL_SAMPLES * 2))

#This lists all population counts/frequencies, allele counts/frequencies, and HWE info
awk -v total_samples="$TOTAL_SAMPLES" -v total_haps="$TOTAL_HAPS" '{print $1, $2, $3, $4, $5, $6, $7, $8, $9, ($9/total_samples), $10, ($10/total_haps), $11, $12, $13, $14, $15, $16, $17, $18, $19}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp10.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp11.txt

#This determines if the SV is in HWE
awk '{
   if ($21 == "NULL") {
       $21 = "HWE_NULL";
   } else if ($21 > 3.84) {
       $21 = "HWE_FALSE";
   } else {
       $21 = "HWE_TRUE";
   }
   print $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, $15, $16, $17, $18, $19, $20, $21
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp11.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp12.txt

##This adds the absolute value of the length of the SV to the start position (DELs are listed as negative lengths in the vcf, but the range is from the start position plus that number of nucleotides)
awk 'BEGIN {OFS="\t"} { if ($7 == "INS") $3 = $2 + 1; else $3 = $2 + ($6 < 0 ? -1*$6 : $6); print }' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp12.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp13.txt

#This annotates each variant with gene info (+/-1kb to try to capture promoters, etc)
bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp13.txt -b $GENES > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp14.txt

#This sorts the file by gene name
sort -k25,25 "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp14.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp15.txt

#This joins the gene name with any OMIM annotations
join -t $'\t' -a 1 -1 25 -2 4 -o '1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,1.10,1.11,1.12,1.13,1.14,1.15,1.16,1.17,1.18,1.19,1.20,1.21,1.25,2.6' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp15.txt "$OMIM_GENES" > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp16.txt

#This condolidates rows that are identical except for the gene and lists the genes together in the last column (ie. a large deletion that spans 50 genes will be consolidated to one row)
awk 'BEGIN{FS=OFS="\t"} {key=$1 FS $2 FS $3 FS $4 FS $5 FS $6 FS $7 FS $8 FS $9 FS $10 FS $11 FS $12 FS $13 FS $14 FS $15 FS $16 FS $17 FS $18 FS $19 FS $20 FS $21} {a[key]=(a[key]!="")?(a[key] ", " $22):(a[key] $22)} {b[key]=(b[key]!="")?(b[key] ", " $23):(b[key] $23)} END {for (k in a) print k, a[k], b[k]}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp16.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp17.txt

#This replaces all empoty fields with a . (important for next command)
awk 'BEGIN {FS = OFS = "\t"} {for (i = 1; i <= NF; i++) if ($i == "") $i = "."} 1' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp17.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp18.txt

#These commands annotates each SV with info (for each bed file), removes unnecessary columns, labels annotated regions as such, and sorts the file (necessary for bedtools intersect)

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp18.txt -b $CODING_REGIONS | \
cut -f1-23,25 | \
awk 'BEGIN {FS=OFS="\t"} $24 == -1 {$24 = "."} $24 != -1 && $24 != "." {$24 = "exonic"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp19.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp19.txt -b $CENTROMERES | \
cut -f1-24,26 | \
awk 'BEGIN {FS=OFS="\t"} $25 == -1 {$25 = "."} $25 != -1 && $25 != "." {$25 = "centromeric"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp20.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp20.txt -b $PERICENTROMERES | \
cut -f1-25,27 | \
awk 'BEGIN {FS=OFS="\t"} $26 == -1 {$26 = "."} $26 != -1 && $26 != "." {$26 = "pericentromeric"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp21.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp21.txt -b $TELOMERES | \
cut -f1-26,28 | \
awk 'BEGIN {FS=OFS="\t"} $27 == -1 {$27 = "."} $27 != -1 && $27 != "." {$27 = "telomeric"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp22.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp22.txt -b $STR | \
cut -f1-27,29 | \
awk 'BEGIN {FS=OFS="\t"} $28 == -1 {$28 = "."} $28 != -1 && $28 != "." {$28 = "STR"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp23.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp23.txt -b $VNTR | \
cut -f1-28,30 | \
awk 'BEGIN {FS=OFS="\t"} $29 == -1 {$29 = "."} $29 != -1 && $29 != "." {$29 = "VNTR"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp24.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp24.txt -b $SEGDUPS | \
cut -f1-29,31 | \
awk 'BEGIN {FS=OFS="\t"} $30 == -1 {$30 = "."} $30 != -1 && $30 != "." {$30 = "segdups"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp25.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp25.txt -b $REPEAT_MASKER | \
cut -f1-30,32 | \
awk 'BEGIN {FS=OFS="\t"} $31 == -1 {$31 = "."} $31 != -1 && $31 != "." {$31 = "repeat"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp26.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp26.txt -b $GAPS | \
cut -f1-31,33 | \
awk 'BEGIN {FS=OFS="\t"} $32 == -1 {$32 = "."} $32 != -1 && $32 != "." {$32 = "gap"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp27.txt

bedtools intersect -wa -wb -loj -a "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp27.txt -b $DEFRABB_HICONF  -f 1.0 | \
cut -f1-32,34 | \
awk 'BEGIN {FS=OFS="\t"} $33 == -1 {$33 = "."} $33 != -1 && $33 != "." {$33 = "hiconf"} 1' | \
sort -u > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp28.txt

# $1-$8 - SV info (Chr, Start_Pos, End_Pos, REF, ALT, SV_Length, SV_Type, 1KGP support)
# $9-$12 - pop/allele freq ($12 = AlleleFreqAll)
# $13-$21 - HWE (GT_homWT, GThet, GT_homVar, p, q, np^2, 2npq, nq^2, HWE)
# $22-$33 - Annotations (gene, OMIM, exonic, cent, peri, tel, str, vntr, segdup, repeat, gaps, hiconf)

#This prints the columns in the correct order
awk 'BEGIN {FS=OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7, $8, $12, $22, $23, $24, $25, $26, $27, $28, $29, $30, $31, $32, $33, $9, $10, $11, $12, $13, $14, $15, $21}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_temp28.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp1.txt

# $1-$8 - SV info (Chr, Start_Pos, End_Pos, REF, ALT, SV_Length, SV_Type, 1KGP support)
# $9 - AlleleFreqAll
# $10-$21 - Annotations (gene, OMIM, exonic, cent, peri, tel, str, vntr, segdup, repeat, gaps, hiconf)
# $22-$25 - pop/allele freq
# $26-$29 - (GT_homWT, GThet, GT_homVar, HWE)

#This fixes strings that are too long for Excel

awk -F'\t' 'BEGIN { OFS="\t" }
{
  if (length($4) > 30000) {
    $4 = "[too long, see vcf]"
  }
  print
}
' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp1.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp2.txt

awk -F'\t' 'BEGIN { OFS="\t" }
{
  if (length($5) > 30000) {
    $5 = "[too long, see vcf]"
  }
  print
}
' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp2.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp3.txt

awk -F'\t' 'BEGIN { OFS="\t" }
{
  if (length($11) > 30000) {
    $11 = "[too long, see vcf]"
  }
  print
}
' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp3.txt > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp4.txt

#This removes "?s" from OMIM annotations (hard to sort)
sed -i 's/?//g' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp4.txt

#This removes empty lines from the bottom of the txt file
sed -i -e :a -e '/^[[:space:]]*$/{$d;N;ba}' "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp4.txt"

#This creates an output file with the correct column headers
echo -e "Chr\tStart_Pos\tEnd_Pos\tREF\tALT\tSV_Length\tSV_Type\tSample_support\tAllele_Freq_ALL\tGenes\tOMIM\tExonic\tCentromeric\tPericentromeric\tTelomeric\tSTR\tVNTR\tSegdup\tRepeat\tGap\tHiConf\tPop_Count_ALL\tPop_Freq_ALL\tAllele_Count_ALL\tAllele_Freq_ALL\tGT_homWT\tGT_het\tGT_homVAR\tHWE" > "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS.txt

#This adds the data to the output file
cat "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp4.txt" >> "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS.txt"

#This makes a vcf of the output
cat <<EOF > "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp.vcf"
##fileformat=VCFv4.2
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr10,length=133797422>
##contig=<ID=chr11,length=135086622>
##contig=<ID=chr12,length=133275309>
##contig=<ID=chr13,length=114364328>
##contig=<ID=chr14,length=107043718>
##contig=<ID=chr15,length=101991189>
##contig=<ID=chr16,length=90338345>
##contig=<ID=chr17,length=83257441>
##contig=<ID=chr18,length=80373285>
##contig=<ID=chr19,length=58617616>
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr21,length=46709983>
##contig=<ID=chr22,length=50818468>
##contig=<ID=chr3,length=198295559>
##contig=<ID=chr4,length=190214555>
##contig=<ID=chr5,length=181538259>
##contig=<ID=chr6,length=170805979>
##contig=<ID=chr7,length=159345973>
##contig=<ID=chr8,length=145138636>
##contig=<ID=chr9,length=138394717>
##contig=<ID=chrM,length=16569>
##contig=<ID=chrX,length=156040895>
##contig=<ID=chrY,length=57227415>
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the SV">
##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="SV length">
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="SV type">
##INFO=<ID=Sample_support,Number=1,Type=String,Description="Samples in the multisample vcf that share the SV">
##INFO=<ID=Allele_Freq_ALL,Number=1,Type=Float,Description="Allele frequency of the SV in the multisample vcf">
##INFO=<ID=Genes,Number=.,Type=String,Description="Genes that the SV overlaps (canonical hg38 coordinates, gencode)">
##INFO=<ID=OMIM,Number=.,Type=String,Description="OMIM phenotypes associated with any gene the SV overlaps (OMIM 8/2023)">
##INFO=<ID=Exonic,Number=1,Type=String,Description="If the SV overlaps exonic regions (canonical hg38 coordinates, gencode)">
##INFO=<ID=Centromeric,Number=1,Type=String,Description="If the SV overlaps centromeric regions (UCSC hg38)">
##INFO=<ID=Pericentromeric,Number=1,Type=String,Description="If the SV intersects with a pericentromeric region (+/-5Mb on either side of UCSC-defined centromere)">
##INFO=<ID=Telomeric,Number=1,Type=String,Description="If the SV intersects with a telomere (5Mb of either end of a chromosome)">
##INFO=<ID=STR,Number=1,Type=String,Description="If the SV intersects with a Short Tandem Repeat region (vamos original motifs, n=148)">
##INFO=<ID=VNTR,Number=1,Type=String,Description="If the SV intersects with a Variable Number Tandem Repeat region (vamos original motifs, n=148)">
##INFO=<ID=Segdup,Number=1,Type=String,Description="If the SV intersects with a segmental duplication (Genome in a Bottle v3.3)">
##INFO=<ID=Repeat,Number=1,Type=String,Description="If the SV intersects with a repeat region (UCSC hg38 repeat masker)">
##INFO=<ID=Gap,Number=1,Type=String,Description="If the SV intersects with a an hg38 gap region (UCSC hg38 mapping and sequencing: gap)">
##INFO=<ID=HiConf,Number=1,Type=String,Description="If the SV is fully contained within a high confidence region (Genome in a Bottle T2TQ100-V1.0_stvar)">
##INFO=<ID=Pop_Count_ALL,Number=1,Type=Integer,Description="Number of individuals in the multisample vcf that carry the SV">
##INFO=<ID=Pop_Freq_ALL,Number=1,Type=Float,Description="Percent of individuals in the multisample vcf that carry the SV">
##INFO=<ID=Allele_Count_ALL,Number=1,Type=Integer,Description="Number of alleles in the cohort with the SV">
##INFO=<ID=Allele_Freq_ALL2,Number=1,Type=Float,Description="(2) Allele frequency of the SV in the multisample vcf">
##INFO=<ID=GT_homWT,Number=1,Type=Integer,Description="Count of homozygous REF genotypes in the multisample vcf">
##INFO=<ID=GT_het,Number=1,Type=Integer,Description="Count of heterozygous genotypes for the SV in the multisample vcf">
##INFO=<ID=GT_homVAR,Number=1,Type=Integer,Description="Count of homozygous ALT genotypes for the SV in the multisample vcf">
##INFO=<ID=HWE,Number=1,Type=String,Description="Hardy-Weinberg Equilibrium result">
EOF

echo -e '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO' >> "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp.vcf"

##Convert updated RESULTS.txt to VCF entries
awk -F"\t" '
BEGIN {
  OFS = "\t"
}
NR>1 {
  chrom = $1
  pos   = $2
  ref   = $4
  alt   = $5

  info = "END=" $3 \
         ";SVLEN=" $6 \
         ";SVTYPE=" $7 \
         ";Sample_support=" $8 \
         ";Allele_Freq_ALL=" $9 \
         ";Genes=" $10 \
         ";OMIM=" $11 \
         ";Exonic=" $12 \
         ";Centromeric=" $13 \
         ";Pericentromeric=" $14 \
         ";Telomeric=" $15 \
         ";STR=" $16 \
         ";VNTR=" $17 \
         ";Segdup=" $18 \
         ";Repeat=" $19 \
         ";Gap=" $20 \
         ";HiConf=" $21 \
         ";Pop_Count_ALL=" $22 \
         ";Pop_Freq_ALL=" $23 \
         ";Allele_Count_ALL=" $24 \
         ";Allele_Freq_ALL2=" $25 \
         ";GT_homWT=" $26 \
         ";GT_het=" $27 \
         ";GT_homVAR=" $28 \
         ";HWE=" $29

  print chrom, pos, ".", ref, alt, ".", ".", info
}' "$UNSOLVED_DIR"/"$QUERY_FILE_NAME"_RESULTS_temp1.txt \
>> "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp.vcf"

bcftools sort "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS_temp.vcf" > "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS.vcf"

bgzip "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS.vcf"
tabix "${UNSOLVED_DIR}/${QUERY_FILE_NAME}_RESULTS.vcf.gz"

##This removes intermediate "temp" files 
find "$UNSOLVED_DIR" -type f -name "*temp*" -exec rm {} \;
