hlakit

#!/usr/bin/env bash

set -uo pipefail

# main

usage() {
    cat << EOF
Usage: $0 --resultdir <dir> --normalbam <bam> --tumorbam <bam> --allelefile <file> --reference <string> --hlakit <path> --tumorpurity <int>

Required arguments:
    -o, --resultdir <dir>         Output directory for results
    -n, --normalbam <bam>         Path to normal BAM file
    -t, --tumorbam <bam>          Path to tumor BAM file
    -a, --allelefile <file>       Path to allele typing result file
    -r, --reference <string>      GRCh reference build used to align the bam files. Accepted values are hg19 and hg38
    -P, --hlakit <path>           Path to hlakit
    -p, --tumorpurity <int>       Tumor purity in fraction

Optional arguements:
    -f, --format <string>         FASTQ format (ILM1.8/ILMFQ/STDFQ; auto-detected if not specified)
    -X, --XMX <int>               Maximum memory for gatk (default: 10)
    -l, --loh <string>            yes/no. yes to run LOH analyis (default: no)
    --WESnormalcoverage <int>     Coverage of normal sample in WES data for running LOH
    --WEStumorcoverage <int>      Coverage of tumor sample in WES data for running LOH
    -@, --threads <int>           Number of threads to use. Use either 4 or 8 (default: 8)
    -s, --shortread <string>      yes/no. Enable short read filter (default: yes)
    -b, --subsettedbam <string>   yes/no. Yes if the input bam files are already subsetted for HLA region  (default: no)
    -h, --help                    Show this help message


Example:
    $0 --resultdir /path/to/resultdirectory --normalbam normalwes.bam --tumorbam tumorwes.bam --allelefile /path/to/allelefile --reference hg19 --hlakit /path/to/hlakit --format ILM1.8 --threads 8 -tumorpurity 0.80 -subsettedbam no
    $0 -o /path/to/resultdirectory -n normal.bam -t tumor.bam -a /path/to/allelefile -r hg19 -P /path/to/hlakit -f ILM1.8 -@ 8 -p 0.80 -b no
EOF
    exit 1
}

# Initialize variables
resultdir=""
normalbam=""
tumorbam=""
allelefile=""
reference=""
loh=no
tumorpurity=""
WESnormalcoverage=""
WEStumorcoverage=""
hlakit=""
format=""
shortread=yes
subsettedbam=no
threads=8
XMX=10

# Parse arguments
while [[ $# -gt 0 ]]; do
    case $1 in
        -o|--resultdir)
            resultdir="$2"
            shift 2
            ;;
        -n|--normalbam)
            normalbam="$2"
            shift 2
            ;;
        -t|--tumorbam)
            tumorbam="$2"
            shift 2
            ;;
        -a|--allelefile)
            allelefile="$2"
            shift 2
            ;;
        -r|--reference)
            reference="$2"
            shift 2
            ;;
        -l|--loh)
            loh="$2"
            shift 2
            ;;
        -p|--tumorpurity)
            tumorpurity="$2"
            shift 2
            ;;
        --WESnormalcoverage)
            WESnormalcoverage="$2"
            shift 2
            ;;
        --WEStumorcoverage)
            WEStumorcoverage="$2"
            shift 2
            ;;
        -X|--XMX)
            XMX="$2"
            shift 2
            ;;
        -P|--hlakit)
            hlakit="$2"
            shift 2
            ;;
        -f|--format)
            format=$2
            shift 2
            ;;
        -@|--threads)
            threads=$2
            shift 2
            ;;
        -s|--shortread)
            shortread=$2
            shift 2
            ;;
        -b|--subsettedbam)
            subsettedbam="$2"
            shift 2
            ;;
        -h|--help)
            usage
            ;;
        *)
            echo "Error: Unknown option $1"
            usage
            ;;
    esac
done

# Validate required arguments
if [ -z "$resultdir" ] || [ -z "$normalbam" ] || [ -z "$tumorbam" ] || [ -z "$allelefile" ] || [ -z "$reference" ] || [ -z "$hlakit" ] || [ -z "$tumorpurity" ]; then
    echo "Error: All arguments are required"
    usage
fi

# Validate reference value
if [ "$reference" != "hg19" ] && [ "$reference" != "hg38" ]; then
    echo "Error: --reference must be 'hg19' or 'hg38' (got: '$reference')."
    exit 1
fi

# Validate shortread and subsettedbam
if [ "$shortread" != "yes" ] && [ "$shortread" != "no" ]; then
    echo "Error: --shortread must be 'yes' or 'no' (got: '$shortread')."
    exit 1
fi

# Validate shortread and subsettedbam
if [ "$subsettedbam" != "yes" ] && [ "$subsettedbam" != "no" ]; then
    echo "Error: --subsettedbam must be 'yes' or 'no' (got: '$subsettedbam')."
    exit 1
fi

# Validate threads and XMX are numeric
if ! [[ "$threads" =~ ^[0-9]+$ ]]; then
    echo "Error: --threads must be a positive integer (got: '$threads')."
    exit 1
fi
if ! [[ "$XMX" =~ ^[0-9]+$ ]]; then
    echo "Error: --XMX must be a positive integer (got: '$XMX')."
    exit 1
fi

# Validate LOH arguments if loh=yes
if [ "$loh" == "yes" ]; then
    if [ -z "$tumorpurity" ]; then
        echo "Error: --tumorpurity is required when --loh yes is set."
        exit 1
    fi
    if [ -z "$WESnormalcoverage" ]; then
        echo "Error: --WESnormalcoverage is required when --loh yes is set."
        exit 1
    fi
    if [ -z "$WEStumorcoverage" ]; then
        echo "Error: --WEStumorcoverage is required when --loh yes is set."
        exit 1
    fi
    if ! [[ "$tumorpurity" =~ ^0(\.[0-9]+)?$|^1(\.0+)?$ ]]; then
        echo "Error: --tumorpurity must be a fraction between 0 and 1 (got: '$tumorpurity'). Do not use a percentage."
        exit 1
    fi
elif [ "$loh" != "no" ]; then
    echo "Error: --loh must be 'yes' or 'no' (got: '$loh')."
    exit 1
fi

# Validate files exist
for file in "$normalbam" "$tumorbam" "$allelefile"; do
    if [ ! -f "$file" ]; then
        echo "Error: File not found: $file"
        exit 1
    elif [ ! -s "$file" ]; then
        echo "Error: File is empty: $file"
        exit 1
    fi
done

# Validate BAM indexes
for bam in "$normalbam" "$tumorbam"; do
    if [ ! -f "${bam}.bai" ] && [ ! -f "${bam%.bam}.bai" ]; then
        echo "Error: BAM index (.bai) not found for $bam — run 'samtools index $bam' before running this script."
        exit 1
    fi
done

# checking if resource files are present
for file in "$hlakit"/resources/kmer.txt "$hlakit"/resources/hla.bed "$hlakit"/resources/hla.dict "$hlakit"/resources/hla.fasta "$hlakit"/resources/hla.fasta.fai "$hlakit"/resources/hla.gtf "$hlakit"/resources/hla.nix "$hlakit"/resources/hlasamheader.sam "$hlakit"/resources/hs_metrics_interval_list.txt
do
    if [ ! -f "$file" ]; then
        echo "Error: Resource file not found: $file"
        exit 1
    elif [ ! -s "$file" ]; then
        echo "Error: Resource file is empty: $file"
        exit 1
    fi
done

# checking if binary files are present
for file in "$hlakit"/binaries/novoalign "$hlakit"/binaries/fastp
do
    if [ ! -f "$file" ]; then
        echo "Error: Binary not found for: $file"
        exit 1
    elif [ ! -s "$file" ]; then
        echo "Error: Binary is empty for: $file"
        exit 1
    fi
done

#checking if software dependencies are installed and in PATH
REQUIRED_TOOLS=(
    "samtools"
    "gatk"
    "Rscript"
)
check_dependencies() {
    local missing=()
    for tool in "${REQUIRED_TOOLS[@]}"; do
        if ! command -v "$tool" &>/dev/null; then
            missing+=("$tool")
        fi
    done
    if [[ ${#missing[@]} -gt 0 ]]; then
        echo "ERROR: Missing required tools:" >&2
        printf '  - %s\n' "${missing[@]}" >&2
        exit 1
    fi
    echo "All required tools found."
}
check_dependencies

# checking if all required R packages are installed
REQUIRED_R_PACKAGES=(
    "Rsamtools"
    "stringr"
    "readr"
    "IRanges"
    "tidyverse"
    "Biostrings"
    "GenomicRanges"
    "optparse"
    "tidyr"
    "dplyr"
    "ggplot2"
)
check_r_packages() {
    local pkg_array
    pkg_array=$(printf '"%s",' "${REQUIRED_R_PACKAGES[@]}")
    pkg_array="c(${pkg_array%,})"  # trim trailing comma

    missing=$(Rscript -e "
        pkgs <- ${pkg_array}
        missing <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
        if (length(missing)) cat(missing, sep = '\n')
    ")

    if [[ -n "$missing" ]]; then
        echo "ERROR: Missing required R packages:" >&2
        sed 's/^/  - /' <<< "$missing" >&2
        exit 1
    fi

    echo "All required R packages found."
}
check_r_packages


# kmer and hla sam header file paths
kmer=$hlakit/resources/kmer.txt
hlasamheader=$hlakit/resources/hlasamheader.sam


# Create output directory
if [ ! -d $resultdir ]; then
    mkdir -p "$resultdir"/result
fi
mkdir -p "$resultdir"/result
resultdir="$resultdir"/result

log_file="${resultdir}/hlakit_$(date +%Y%m%d_%H%M%S).log"
mkfifo /tmp/logpipe_$$
tee -a "$log_file" < /tmp/logpipe_$$ &
exec > /tmp/logpipe_$$ 2>&1

# color-coded info, warning, and error
RED='\033[0;31m'
YELLOW='\033[1;33m'
GREEN='\033[0;32m'
NC='\033[0m'

log_info() {
    local msg="[$(date +'%Y-%m-%d %H:%M:%S')] [INFO] $*"
    echo -e "${GREEN}${msg}${NC}" >&2
    # echo "$msg" >> "$log_file"
}

log_warn() {
    local msg="[$(date +'%Y-%m-%d %H:%M:%S')] [WARNING] $*"
    echo -e "${YELLOW}${msg}${NC}" >&2
    # echo "$msg" >> "$log_file"
}

log_error() {
    local msg="[$(date +'%Y-%m-%d %H:%M:%S')] [ERROR] $*"
    echo -e "${RED}${msg}${NC}" >&2
    # echo "$msg" >> "$log_file"
}


run_command() {
    local cmd="$*"
    log_info "Running: $cmd"
    
    # Only log stderr, not stdout
    if eval "$cmd" 2>> "$log_file"; then
        log_info "Finished: $cmd"
    else
        log_error "Command failed: $cmd"
        return 1
    fi
}

log_info "Logging session in: $log_file"

log_info "####################################################"
log_info "################## BEGIN ANALYSIS ##################"
log_info "####################################################"

echo >> $allelefile
sort $allelefile | uniq | sed -e '1{/^$/d}' -e '$a\\' | tr -d '\r' > $resultdir/allelefile.tmp.txt && mv $resultdir/allelefile.tmp.txt $allelefile
log_info "Alleles in input allelefile:"
run_command cat $allelefile

:>"$resultdir/zygosity.txt"
for gene in hla_a hla_b hla_c; do
    count=$(grep -c "^${gene}" "$allelefile" || true)
    if [ "$count" -eq 2 ]; then
        status="Heterozygous"
    elif [ "$count" -eq 1 ]; then
        status="Homozygous"
    else
        echo "WARNING: unexpected allele count for ${gene}: ${count}" >&2
        status="Unknown"
    fi
    printf "%s\t%s\n" "$gene" "$status"
done >> "$resultdir/zygosity.txt"

log_info "####################################################"
log_info "################# BAM SUBSETTING ###################"
log_info "####################################################"


log_info "SUBSETTING NORMAL BAM"

if [[ "$subsettedbam" == "yes" ]]; then
    log_info "Converting HLA region reads from Normal SAM to FASTQ ..."
    run_command samtools sort -o ${normalbam/.bam/.sorted.bam} $normalbam
    run_command samtools index ${normalbam/.bam/.sorted.bam}
    run_command gatk SamToFastq -I ${normalbam/.bam/.sorted.bam} -F $resultdir/normal_subsetted.1.fastq -F2 $resultdir/normal_subsetted.2.fastq --VALIDATION_STRINGENCY SILENT

    log_info "Removing unpaired reads from Normal chr6 HLA region Normal FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/normal_subsetted.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/normal_subsetted.2.fastq

    if [[ ! -s "$resultdir/normal_subsetted.1.fastq" || ! -s "$resultdir/normal_subsetted.2.fastq" ]]; then
        log_error "Merged Normal FASTQ files are missing or empty after merge_fastq.sh: $resultdir/normal_subsetted.1.fastq / .2.fastq — check input Normal bam file."
        exit 1
    fi

    log_info "Converting HLA region reads from Tumor SAM to FASTQ ..."
    run_command samtools sort -o ${tumorbam/.bam/.sorted.bam} $tumorbam
    run_command samtools index ${tumorbam/.bam/.sorted.bam}
    run_command gatk SamToFastq -I ${tumorbam/.bam/.sorted.bam} -F $resultdir/tumor_subsetted.1.fastq -F2 $resultdir/tumor_subsetted.2.fastq --VALIDATION_STRINGENCY SILENT

    log_info "Removing unpaired reads from tumor chr6 HLA region tumor FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/tumor_subsetted.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/tumor_subsetted.2.fastq

    if [[ ! -s "$resultdir/tumor_subsetted.1.fastq" || ! -s "$resultdir/tumor_subsetted.2.fastq" ]]; then
        log_error "Merged tumor FASTQ files are missing or empty after merge_fastq.sh: $resultdir/tumor_subsetted.1.fastq / .2.fastq — check input Tumor bam file."
        exit 1
    fi
else
    log_info "Extracting k-mer reads from Normal bam ($normalbam) ..."
    run_command samtools view -H $normalbam > $resultdir/normal_kmer.sam
    empty_size=$(stat -c%s $resultdir/normal_kmer.sam)
    run_command samtools view --threads $threads $normalbam | grep -F -f $kmer >> $resultdir/normal_kmer.sam
    run_command samtools view -bS --threads $threads -o $resultdir/normal_kmer.bam $resultdir/normal_kmer.sam
    if [[ $? -ne 0 || $(stat -c%s $resultdir/normal_kmer.sam) -le $empty_size ]]; then
        log_warn "0 reads mapped to k-mer in the Normal bam ($normalbam)."
    fi

    log_info "Converting k-mer reads from Normal SAM to FASTQ ..."
    run_command gatk SamToFastq -I $resultdir/normal_kmer.bam -F $resultdir/normal_kmer.1.fastq -F2 $resultdir/normal_kmer.2.fastq --VALIDATION_STRINGENCY SILENT
    if [[ $? -ne 0 || $(stat -c%s $resultdir/normal_kmer.1.fastq) -le 1 || $(stat -c%s $resultdir/normal_kmer.2.fastq) -le 1 ]]; then
        log_warn "Normal $resultdir/normal_kmer.1.fastq or $resultdir/normal_kmer.2.fastq file is empty."
    fi

    log_info "Removing unpaired reads from Normal k-mer FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/normal_kmer.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/normal_kmer.2.fastq
    if [[ $? -ne 0 || $(stat -c%s $resultdir/normal_kmer.1.fastq) -le 1 || $(stat -c%s $resultdir/normal_kmer.2.fastq) -le 1 ]]; then
        log_warn "Normal $resultdir/normal_kmer.1.fastq or $resultdir/normal_kmer.2.fastq file is empty."
    fi

    log_info "Extracting HLA regions from chr6 of $normalbam ..."
    run_command samtools view --threads $threads -H $normalbam > $resultdir/chr6region.normal.sam
    empty_size=$(stat -c%s $resultdir/chr6region.normal.sam)

    chr_prefix=$(samtools view -H $normalbam | awk '/^@SQ/ && /SN:chr6/{found=1} END{print (found ? "chr" : "")}')
    if [ $reference == "hg38" ]; then
        regions="${chr_prefix}6:29941260-29945884 ${chr_prefix}6:31353872-31357187 ${chr_prefix}6:31268749-31272105"
    elif [ $reference == "hg19" ]; then
        regions="${chr_prefix}6:29909037-29913661 ${chr_prefix}6:31321649-31324964 ${chr_prefix}6:31236526-31239869"
    fi
    run_command samtools view --threads $threads $normalbam $regions >> $resultdir/chr6region.normal.sam

    if [[ $? -ne 0 || $(stat -c%s $resultdir/chr6region.normal.sam) -le "$empty_size" ]]; then
        log_warn "0 reads mapped to chr6 HLA regions in the Normal bam ($normalbam)."
    fi

    log_info "Creating subsetted Normal BAM with kmer and chr6 region reads ..."
    # Merge SAMs
    cat $resultdir/normal_kmer.sam <(samtools view $resultdir/chr6region.normal.sam) | \
        awk '(/^@/ && !seen_header[$0]++) || (!/^@/ && !seen[$1":"$2]++)' \
        > $resultdir/normal_subsetted.${reference}.sam \
        || log_warn "Failed to merge SAM files"

    [ -s $resultdir/normal_subsetted.${reference}.sam ] || log_warn "Merged SAM is empty"
    # SAM to BAM
    run_command samtools view -bS -o $resultdir/normal_subsetted.${reference}.bam $resultdir/normal_subsetted.${reference}.sam
    rm -f $resultdir/normal_subsetted.${reference}.sam
    # Coordinate sort
    run_command gatk SortSam -I $resultdir/normal_subsetted.${reference}.bam -O $resultdir/normal_subsetted.${reference}.coordsort.bam --SORT_ORDER coordinate
    rm -f $resultdir/normal_subsetted.${reference}.bam
    # Mark duplicates
    run_command gatk --java-options -Xmx${XMX}g MarkDuplicates -I $resultdir/normal_subsetted.${reference}.coordsort.bam -O $resultdir/normal_subsetted.${reference}.dedup.bam -M $resultdir/normal_subsetted.${reference}.dedup.metrics.txt --REMOVE_DUPLICATES true
    # log_info "MarkDuplicates metrics:"; echo $resultdir/normal_subsetted.${reference}.dedup.metrics.txt
    rm -f $resultdir/normal_subsetted.${reference}.coordsort.bam $resultdir/normal_subsetted.${reference}.dedup.metrics.txt
    run_command samtools sort -o $resultdir/normal_subsetted.${reference}.bam $resultdir/normal_subsetted.${reference}.dedup.bam
    rm -f $resultdir/normal_subsetted.${reference}.dedup.bam
    run_command samtools index $resultdir/normal_subsetted.${reference}.bam
    [ -f $resultdir/normal_subsetted.${reference}.bam.bai ] || log_warn "BAM index (.bai) not created"


    log_info "Converting chr6 HLA region reads from Normal SAM to FASTQ ..."
    run_command samtools view -bS -o $resultdir/chr6region.normal.bam $resultdir/chr6region.normal.sam
    run_command gatk SamToFastq -I $resultdir/chr6region.normal.bam -F $resultdir/chr6region.normal.1.fastq -F2 $resultdir/chr6region.normal.2.fastq --VALIDATION_STRINGENCY SILENT

    log_info "Removing unpaired reads from Normal chr6 HLA region Normal FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/chr6region.normal.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/chr6region.normal.2.fastq

    log_info "Merging reads from Normal k-mer and chr6 HLA region ..."
    run_command bash "$hlakit"/merge_fastq.sh -k $resultdir/normal_kmer -c $resultdir/chr6region.normal -o $resultdir/normal_subsetted

    if [[ ! -s "$resultdir/normal_subsetted.1.fastq" || ! -s "$resultdir/normal_subsetted.2.fastq" ]]; then
        log_error "Merged Normal FASTQ files are missing or empty after merge_fastq.sh: $resultdir/normal_subsetted.1.fastq / .2.fastq — check k-mer and chr6 extraction steps above."
        exit 1
    fi


    log_info "SUBSETTING TUMOR BAM"
    log_info "Extracting k-mer reads from Tumor bam ($tumorbam) ..."
    run_command samtools view -H $tumorbam > $resultdir/tumor_kmer.sam
    empty_size=$(stat -c%s $resultdir/tumor_kmer.sam)
    run_command samtools view --threads $threads $tumorbam | grep -F -f $kmer >> $resultdir/tumor_kmer.sam
    run_command samtools view -bS --threads $threads -o $resultdir/tumor_kmer.bam $resultdir/tumor_kmer.sam
    if [[ $? -ne 0 || $(stat -c%s $resultdir/tumor_kmer.sam) -le $empty_size ]]; then
        log_warn "0 reads mapped to k-mer in the Tumor bam ($tumorbam)."
    fi


    log_info "Converting k-mer reads from Tumor SAM to FASTQ ..."
    run_command gatk SamToFastq -I $resultdir/tumor_kmer.bam -F $resultdir/tumor_kmer.1.fastq -F2 $resultdir/tumor_kmer.2.fastq --VALIDATION_STRINGENCY SILENT
    if [[ $? -ne 0 || $(stat -c%s $resultdir/tumor_kmer.1.fastq) -le 1 || $(stat -c%s $resultdir/tumor_kmer.2.fastq) -le 1 ]]; then
        log_warn "Tumor $resultdir/tumor_kmer.1.fastq or $resultdir/tumor_kmer.2.fastq file is empty."
    fi

    log_info "Removing unpaired reads from Tumor k-mer FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/tumor_kmer.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/tumor_kmer.2.fastq
    if [[ $? -ne 0 || $(stat -c%s $resultdir/tumor_kmer.1.fastq) -le 1 || $(stat -c%s $resultdir/tumor_kmer.2.fastq) -le 1 ]]; then
        log_warn "Tumor $resultdir/tumor_kmer.1.fastq or $resultdir/tumor_kmer.2.fastq file is empty."
    fi

    log_info "Extracting HLA regions from $tumorbam ..."
    run_command samtools view --threads $threads -H $tumorbam > $resultdir/chr6region.tumor.sam
    empty_size=$(stat -c%s $resultdir/chr6region.tumor.sam)


    chr_prefix=$(samtools view -H $tumorbam | awk '/^@SQ/ && /SN:chr6/{found=1} END{print (found ? "chr" : "")}')
    if [ $reference == "hg38" ]; then
        regions="${chr_prefix}6:29941260-29945884 ${chr_prefix}6:31353872-31357187 ${chr_prefix}6:31268749-31272105"
    elif [ $reference == "hg19" ]; then
        regions="${chr_prefix}6:29909037-29913661 ${chr_prefix}6:31321649-31324964 ${chr_prefix}6:31236526-31239869"
    fi
    run_command samtools view --threads $threads $tumorbam $regions >> $resultdir/chr6region.tumor.sam

    if [[ $? -ne 0 || $(stat -c%s $resultdir/chr6region.tumor.sam) -le "$empty_size" ]]; then
        log_warn "0 reads mapped to chr6 HLA Class I regions in the Tumor bam ($tumorbam)."
    fi


    log_info "Creating subsetted Tumor BAM with kmer and chr6 region reads ..."
    # Merge SAMs
    cat $resultdir/tumor_kmer.sam <(samtools view $resultdir/chr6region.tumor.sam) | \
        awk '(/^@/ && !seen_header[$0]++) || (!/^@/ && !seen[$1":"$2]++)' \
        > $resultdir/tumor_subsetted.${reference}.sam \
        || log_warn "Failed to merge SAM files"

    [ -s $resultdir/tumor_subsetted.${reference}.sam ] || log_warn "Merged SAM is empty"
    # SAM to BAM
    run_command samtools view -bS -o $resultdir/tumor_subsetted.${reference}.bam $resultdir/tumor_subsetted.${reference}.sam
    rm -f $resultdir/tumor_subsetted.${reference}.sam
    # Coordinate sort
    run_command gatk SortSam -I $resultdir/tumor_subsetted.${reference}.bam -O $resultdir/tumor_subsetted.${reference}.coordsort.bam --SORT_ORDER coordinate
    rm -f $resultdir/tumor_subsetted.${reference}.bam
    # Mark duplicates
    run_command gatk --java-options -Xmx${XMX}g MarkDuplicates -I $resultdir/tumor_subsetted.${reference}.coordsort.bam -O $resultdir/tumor_subsetted.${reference}.dedup.bam -M $resultdir/tumor_subsetted.${reference}.dedup.metrics.txt --REMOVE_DUPLICATES true
    # log_info "MarkDuplicates metrics:"; echo $resultdir/tumor_subsetted.${reference}.dedup.metrics.txt
    rm -f $resultdir/tumor_subsetted.${reference}.coordsort.bam $resultdir/tumor_subsetted.${reference}.dedup.metrics.txt
    run_command samtools sort -o $resultdir/tumor_subsetted.${reference}.bam $resultdir/tumor_subsetted.${reference}.dedup.bam
    rm -f $resultdir/tumor_subsetted.${reference}.dedup.bam
    run_command samtools index $resultdir/tumor_subsetted.${reference}.bam
    [ -f $resultdir/tumor_subsetted.${reference}.bam.bai ] || log_warn "BAM index (.bai) not created"


    log_info "Converting chr6 HLA region reads from Tumor SAM to FASTQ ..."
    run_command samtools view -bS -o $resultdir/chr6region.tumor.bam $resultdir/chr6region.tumor.sam
    run_command gatk SamToFastq -I $resultdir/chr6region.tumor.bam -F $resultdir/chr6region.tumor.1.fastq -F2 $resultdir/chr6region.tumor.2.fastq --VALIDATION_STRINGENCY SILENT

    log_info "Removing unpaired reads from Tumor chr6 HLA region Tumor FASTQ ..."
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/chr6region.tumor.1.fastq
    run_command bash "$hlakit"/paired_fastq.sh -f $resultdir/chr6region.tumor.2.fastq

    log_info "Merging reads from Tumor k-mer and chr6 HLA region ..."
    run_command bash "$hlakit"/merge_fastq.sh -k $resultdir/tumor_kmer -c $resultdir/chr6region.tumor -o $resultdir/tumor_subsetted

    if [[ ! -s "$resultdir/tumor_subsetted.1.fastq" || ! -s "$resultdir/tumor_subsetted.2.fastq" ]]; then
        log_error "Merged Tumor FASTQ files are missing or empty after merge_fastq.sh: $resultdir/tumor_subsetted.1.fastq / .2.fastq — check k-mer and chr6 extraction steps above."
        exit 1
    fi
fi


log_info "####################################################"
log_info "#################### ALIGNMENT #####################"
log_info "####################################################"


log_info "Cleaning Subsetted Normal fastq files using fastp ..."
run_command "$hlakit"/binaries/fastp -i $resultdir/normal_subsetted.1.fastq -I $resultdir/normal_subsetted.2.fastq -o $resultdir/clean.normal_subsetted.1.fastq -O $resultdir/clean.normal_subsetted.2.fastq -j $resultdir/normal_subsetted.json -h $resultdir/normal_subsetted.html
mv $resultdir/clean.normal_subsetted.1.fastq $resultdir/normal_subsetted.1.fastq
mv $resultdir/clean.normal_subsetted.2.fastq $resultdir/normal_subsetted.2.fastq

if [[ ! -s "$resultdir/normal_subsetted.1.fastq" || ! -s "$resultdir/normal_subsetted.2.fastq" ]]; then
    log_error "Normal FASTQ files are empty after fastp cleaning. Check fastp log: $resultdir/normal_subsetted.json"
    exit 1
fi

log_info "Aligning Normal reads to HLA library ..."
run_command cat "$hlasamheader" > $resultdir/novoalign_normal.sam
if [ ! -z $format ]; then 
    run_command bash "$hlakit"/runnovoalign.sh --f1 $resultdir/normal_subsetted.1.fastq --f2 $resultdir/normal_subsetted.2.fastq --threads $threads --nixFile $hlakit/resources/hla.nix --output $resultdir/novoalign_normal_reads.sam --hlakit $hlakit --format $format
else 
    run_command bash "$hlakit"/runnovoalign.sh --f1 $resultdir/normal_subsetted.1.fastq --f2 $resultdir/normal_subsetted.2.fastq --threads $threads --nixFile $hlakit/resources/hla.nix --output $resultdir/novoalign_normal_reads.sam --hlakit $hlakit
fi

if [[ ! -s "$resultdir/novoalign_normal_reads.sam" ]]; then
    log_error "Normal novoalign output SAM is missing or empty: $resultdir/novoalign_normal_reads.sam — novoalign may have failed or produced 0 alignments."
    exit 1
fi

run_command cat $resultdir/novoalign_normal_reads.sam >> $resultdir/novoalign_normal.sam
run_command samtools view -bS -o $resultdir/novoalign_normal.bam $resultdir/novoalign_normal.sam

log_info "Deduplicating HLA-aligned Normal bam ..."
run_command gatk SortSam -I $resultdir/novoalign_normal.bam -O $resultdir/novoalign_normal.coordsort.bam --SORT_ORDER coordinate
normalsamplename=$(basename $normalbam)
run_command gatk AddOrReplaceReadGroups -I $resultdir/novoalign_normal.coordsort.bam -O $resultdir/novoalign_normal.coordsort.RG.bam -ID "$normalsamplename".normal -LB library1 -PL illumina -PU unit1 -SM "$normalsamplename".normal
run_command gatk --java-options -Xmx${XMX}g MarkDuplicates -I $resultdir/novoalign_normal.coordsort.RG.bam -O $resultdir/novoalign_normal.coordsort.dedup.RG.bam -M $resultdir/novoalign_normal.coordsort.dedup.RG.metrics.txt --REMOVE_DUPLICATES true

if [[ ! -s "$resultdir/novoalign_normal.coordsort.dedup.RG.bam" ]]; then
    log_error "Normal deduplicated BAM is missing or empty: $resultdir/novoalign_normal.coordsort.dedup.RG.bam — MarkDuplicates may have failed or removed all reads."
    exit 1
fi

run_command samtools index $resultdir/novoalign_normal.coordsort.dedup.RG.bam

log_info "Cleaning Subsetted Tumor fastq files using fastp ..."
run_command "$hlakit"/binaries/fastp -i $resultdir/tumor_subsetted.1.fastq -I $resultdir/tumor_subsetted.2.fastq -o $resultdir/clean.tumor_subsetted.1.fastq -O $resultdir/clean.tumor_subsetted.2.fastq -j $resultdir/tumor_subsetted.json -h $resultdir/tumor_subsetted.html
mv $resultdir/clean.tumor_subsetted.1.fastq $resultdir/tumor_subsetted.1.fastq
mv $resultdir/clean.tumor_subsetted.2.fastq $resultdir/tumor_subsetted.2.fastq

if [[ ! -s "$resultdir/tumor_subsetted.1.fastq" || ! -s "$resultdir/tumor_subsetted.2.fastq" ]]; then
    log_error "Tumor FASTQ files are empty after fastp cleaning. Check fastp log: $resultdir/tumor_subsetted.json"
    exit 1
fi

log_info "Aligning Tumor reads to HLA library ..."
run_command cat "$hlasamheader" > $resultdir/novoalign_tumor.sam
if [ ! -z $format ]; then
    run_command bash "$hlakit"/runnovoalign.sh --f1 $resultdir/tumor_subsetted.1.fastq --f2 $resultdir/tumor_subsetted.2.fastq --threads $threads --nixFile $hlakit/resources/hla.nix --output $resultdir/novoalign_tumor_reads.sam --hlakit $hlakit --format $format
else 
    run_command bash "$hlakit"/runnovoalign.sh --f1 $resultdir/tumor_subsetted.1.fastq --f2 $resultdir/tumor_subsetted.2.fastq --threads $threads --nixFile $hlakit/resources/hla.nix --output $resultdir/novoalign_tumor_reads.sam --hlakit $hlakit
fi

if [[ ! -s "$resultdir/novoalign_tumor_reads.sam" ]]; then
    log_error "Tumor novoalign output SAM is missing or empty: $resultdir/novoalign_tumor_reads.sam — novoalign may have failed or produced 0 alignments."
    exit 1
fi

run_command cat $resultdir/novoalign_tumor_reads.sam >> $resultdir/novoalign_tumor.sam
run_command samtools view -bS -o $resultdir/novoalign_tumor.bam $resultdir/novoalign_tumor.sam

log_info "Deduplicating HLA-aligned Tumor bam ..."
run_command gatk SortSam -I $resultdir/novoalign_tumor.bam -O $resultdir/novoalign_tumor.coordsort.bam --SORT_ORDER coordinate
tumorsamplename=$(basename "$tumorbam")
run_command gatk AddOrReplaceReadGroups -I $resultdir/novoalign_tumor.coordsort.bam -O $resultdir/novoalign_tumor.coordsort.RG.bam -ID "$tumorsamplename".tumor -LB library1 -PL illumina -PU unit1 -SM "$tumorsamplename".tumor
run_command gatk --java-options -Xmx${XMX}g MarkDuplicates -I $resultdir/novoalign_tumor.coordsort.RG.bam -O $resultdir/novoalign_tumor.coordsort.dedup.RG.bam -M $resultdir/novoalign_tumor.coordsort.dedup.RG.metrics.txt --REMOVE_DUPLICATES true

if [[ ! -s "$resultdir/novoalign_tumor.coordsort.dedup.RG.bam" ]]; then
    log_error "Tumor deduplicated BAM is missing or empty: $resultdir/novoalign_tumor.coordsort.dedup.RG.bam — MarkDuplicates may have failed or removed all reads."
    exit 1
fi

run_command samtools index $resultdir/novoalign_tumor.coordsort.dedup.RG.bam


log_info "####################################################"
log_info "################### MUTATION CALLING ###############"
log_info "####################################################"

log_info "Verifying if the typed HLA alleles have genomic sequence in the HLA library. If genomic sequence is not available for an allele, hlakit will try to find an allele at next field level. If no match is found, hlakit will skip that allele ..."
var=`cat $allelefile`

alleles=($var)
updated_alleles=()

for allele in "${alleles[@]}"; do
    log_info "Working on allele: $allele ..."
    if [ ! -z "$(cat "$hlasamheader" | grep -P "$allele\t")" ]
    then
        updated_alleles+=("$allele")
    elif [ -z "$(cat "$hlasamheader" | grep -P "$allele\t")" ]
    then
        log_info "Exact match for allele: $allele not found in the HLA reference. Trying to find a match at 6 or 8 field level ..."
        if [ $(cat "$hlasamheader" | grep -P "$allele" | wc -l ) -gt 0 ]
        then
            allele=$(cat "$hlasamheader" | grep -P "$allele"_ | head -1 | awk '{print $2}' FS="\t" | awk '{print $2}' FS=":")
            updated_alleles+=("$allele")
        elif [ $(cat "$hlasamheader" | grep -P "$allele"_ | wc -l ) -gt 0 ]
        then
            allele=$(cat "$hlasamheader" | grep -P "$allele"_ | head -1 | awk '{print $2}' FS="\t" | awk '{print $2}' FS=":")
            updated_alleles+=("$allele")
        else
            log_info "Match not found at 6 or 8 field level. hlakit will ignore $allele ..."
            updated_alleles+=("")
        fi
    fi
done

var="${updated_alleles[*]}"
echo "$var" | tr ' ' '\n' > $resultdir/allelelist_updated.txt

log_info "Updated alleles: $var"

# Check that at least one allele resolved successfully
resolved_alleles=()
for a in "${updated_alleles[@]}"; do
    if [ -n "$a" ]; then resolved_alleles+=("$a"); fi
done
if [ ${#resolved_alleles[@]} -eq 0 ]; then
    log_error "No alleles could be resolved against the HLA reference. Check that the allele file contains valid HLA allele names and that the HLA SAM header is correct."
    exit 1
fi
if [ ${#resolved_alleles[@]} -lt ${#alleles[@]} ]; then
    log_warn "$((${#alleles[@]} - ${#resolved_alleles[@]})) allele(s) out of ${#alleles[@]} could not be resolved and will be skipped."
fi

log_info "Subsetting bam aligned to HLA library for the typed alleles ..."
for BAM in novoalign_normal.coordsort.dedup.RG.bam novoalign_tumor.coordsort.dedup.RG.bam
do
    log_info "Subsetting $BAM ..."
    # run_command samtools view -H $resultdir/"$BAM" > $resultdir/"${BAM/.bam/.typedsubset.sam}"
    run_command samtools view $resultdir/"$BAM" ${var[@]} > $resultdir/"${BAM/.bam/.typedsubset.txt}"

    if [ ! -s "$resultdir/${BAM/.bam/.typedsubset.txt}" ]; then
        log_warn "Typed subset file is empty for $BAM: $resultdir/${BAM/.bam/.typedsubset.txt} — no reads aligned to any of the resolved alleles."
    fi

    log_info "Cleaning alignments for SNP+INDEL detection bam files..."
    run_command bash $hlakit/filterReads.sh -i $resultdir/"${BAM/.bam/.typedsubset.txt}" -o $resultdir/"${BAM/.bam/.typedsubset.cleanSNP.txt}" -m 1

    log_info "Cleaning alignments for DNP detection bam files ..."
    run_command bash $hlakit/filterReads.sh -i $resultdir/"${BAM/.bam/.typedsubset.txt}" -o $resultdir/"${BAM/.bam/.typedsubset.cleanDNP.txt}" -m 2

    if [ ! -s "$resultdir/${BAM/.bam/.typedsubset.cleanSNP.txt}" ]; then
        log_warn "Clean SNP alignment file is empty for $BAM — filterReads.sh may have filtered all reads."
    fi
    if [ ! -s "$resultdir/${BAM/.bam/.typedsubset.cleanDNP.txt}" ]; then
        log_warn "Clean DNP alignment file is empty for $BAM — filterReads.sh may have filtered all reads."
    fi

    #Remove unmapped reads if any. They cause errors in hs metrics
    for file in $resultdir/"${BAM/.bam/.typedsubset.cleanSNP.txt}" $resultdir/"${BAM/.bam/.typedsubset.cleanDNP.txt}"
    do
        faulty_reads=$(awk '{print $1"\t"$6}' $file | egrep "\*" | awk '{print $1}')  || true
        if [ -n "$faulty_reads" ]
        then
            egrep -v "$faulty_reads" $file > ${file/.txt/.temp.txt} && mv ${file/.txt/.temp.txt} $file
        fi
    done
done


normalSNP=novoalign_normal.coordsort.dedup.RG.typedsubset.cleanSNP.txt
tumorSNP=novoalign_tumor.coordsort.dedup.RG.typedsubset.cleanSNP.txt
normalDNP=novoalign_normal.coordsort.dedup.RG.typedsubset.cleanDNP.txt
tumorDNP=novoalign_tumor.coordsort.dedup.RG.typedsubset.cleanDNP.txt


for muttype in SNP DNP
do
    for allele in $var
    do
        # Skip empty allele slots (unresolved alleles)
        if [ -z "$allele" ]; then
            log_warn "Skipping empty allele entry for muttype $muttype."
            continue
        fi

        normal="normal${muttype}" ; normal="${!normal}"
        tumor="tumor${muttype}" ; tumor="${!tumor}"
        
        normalBam=${normal/.txt/.bam}
        tumorBam=${tumor/.txt/.bam}

        log_info "Subsetting tumor and normal bam files to find ${muttype}s in allele $allele ..."

        run_command head -2 "$hlasamheader" > $resultdir/${allele}.normal.${muttype}.sam
        run_command cat "$hlasamheader" | grep -P "$allele\t" >> $resultdir/${allele}.normal.${muttype}.sam
        run_command cp $resultdir/${allele}.normal.${muttype}.sam $resultdir/${allele}.tumor.${muttype}.sam

        printf "@RG\tID:"$normalsamplename".normal\tPL:illumina\tPU:unit1\tLB:library1\tSM:"$normalsamplename".normal\n" >> $resultdir/${allele}.normal.${muttype}.sam
        run_command cat $resultdir/$normal | grep -P "\t$allele\t" | sort -k1 > $resultdir/${allele}.normal.${muttype}.txt

        if [ ! -s $resultdir/${allele}.normal.${muttype}.txt ]; then
            log_warn "$resultdir/${allele}.normal.${muttype}.txt is empty — no normal reads found for allele $allele ($muttype). Mutect2 will likely fail for this allele."
        fi

        run_command bash $hlakit/paired_reads_bam.sh -i $resultdir/${allele}.normal.${muttype}.txt
        run_command cat $resultdir/${allele}.normal.${muttype}.txt >> $resultdir/${allele}.normal.${muttype}.sam
        run_command samtools view -bS -o $resultdir/${allele}.normal.${muttype}.bam $resultdir/${allele}.normal.${muttype}.sam
        run_command samtools sort -o $resultdir/${allele}.normal.${muttype}.sorted.bam $resultdir/${allele}.normal.${muttype}.bam
        run_command samtools index $resultdir/${allele}.normal.${muttype}.sorted.bam

        printf "@RG\tID:"$tumorsamplename".tumor\tPL:illumina\tPU:unit1\tLB:library1\tSM:"$tumorsamplename".tumor\n" >> $resultdir/${allele}.tumor.${muttype}.sam
        run_command cat $resultdir/$tumor | grep -P "\t$allele\t" | sort -k1 > $resultdir/${allele}.tumor.${muttype}.txt

        if [ ! -s $resultdir/${allele}.tumor.${muttype}.txt ]; then
            log_warn "$resultdir/${allele}.tumor.${muttype}.txt is empty — no tumor reads found for allele $allele ($muttype). Mutect2 will likely fail for this allele."
        fi

        run_command bash $hlakit/paired_reads_bam.sh -i $resultdir/${allele}.tumor.${muttype}.txt
        run_command cat $resultdir/${allele}.tumor.${muttype}.txt >> $resultdir/${allele}.tumor.${muttype}.sam
        run_command samtools view -bS -o $resultdir/${allele}.tumor.${muttype}.bam $resultdir/${allele}.tumor.${muttype}.sam
        run_command samtools sort -o $resultdir/${allele}.tumor.${muttype}.sorted.bam $resultdir/${allele}.tumor.${muttype}.bam
        run_command samtools index $resultdir/${allele}.tumor.${muttype}.sorted.bam

        run_command gatk AddOrReplaceReadGroups -I $resultdir/${allele}.normal.${muttype}.sorted.bam -O $resultdir/${allele}.normal.${muttype}.sorted.RG.bam -ID "$normalsamplename".normal -LB library1 -PL illumina -PU unit1 -SM "$normalsamplename".normal
        run_command gatk AddOrReplaceReadGroups -I $resultdir/${allele}.tumor.${muttype}.sorted.bam -O $resultdir/${allele}.tumor.${muttype}.sorted.RG.bam -ID "$tumorsamplename".tumor -LB library1 -PL illumina -PU unit1 -SM "$tumorsamplename".tumor

        run_command mv $resultdir/${allele}.normal.${muttype}.sorted.RG.bam $resultdir/${allele}.normal.${muttype}.bam
        run_command mv $resultdir/${allele}.tumor.${muttype}.sorted.RG.bam $resultdir/${allele}.tumor.${muttype}.bam
        run_command samtools index $resultdir/${allele}.normal.${muttype}.bam
        run_command samtools index $resultdir/${allele}.tumor.${muttype}.bam

        log_info "Running Mutect2 on $allele to find $muttype ..."
        run_command gatk Mutect2 \
        -I $resultdir/${allele}.normal.${muttype}.bam \
        -I $resultdir/${allele}.tumor.${muttype}.bam \
        -normal "$normalsamplename".normal \
        -O $resultdir/${allele}.${muttype}.prefilt.vcf \
        -R $hlakit/resources/hla.fasta

        if [ ! -f "$resultdir/${allele}.${muttype}.prefilt.vcf" ]; then
            log_warn "Mutect2 VCF not produced for allele $allele ($muttype): $resultdir/${allele}.${muttype}.prefilt.vcf — Mutect2 may have failed silently."
        fi

        log_info "Running FilterMutectCalls on $allele to filter $muttype ..."
        run_command gatk FilterMutectCalls \
        -V $resultdir/${allele}.${muttype}.prefilt.vcf \
        -O $resultdir/${allele}.${muttype}.vcf \
        -R $hlakit/resources/hla.fasta

    done
done


log_info "Running MAPQ correction ..."
run_command bash $hlakit/mapq_correction.sh -o $resultdir -a $resultdir/allelelist_updated.txt -t $threads

log_info "Calculating coverage of each allele in normal and tumor bam files ..."
run_command bash $hlakit/calculate_coverage.sh -o $resultdir -a $resultdir/allelelist_updated.txt -p $hlakit
printf "allele\tnormal\ttumor\n" > $resultdir/coverage.txt
for allele in $var
do
    if [ -z "$allele" ]; then continue; fi

    if [ ! -f "$resultdir/${allele}.coverage.txt" ]; then
        log_warn "Coverage file not found for allele $allele: $resultdir/${allele}.coverage.txt — coverage for this allele will be missing from $resultdir/coverage.txt."
        continue
    fi

    normal=$(awk 'NR==1{print $2}' $resultdir/${allele}.coverage.txt)
    tumor=$(awk 'NR==2{print $2}' $resultdir/${allele}.coverage.txt)

    if [ -z "$normal" ] || [ -z "$tumor" ]; then
        log_warn "Could not parse normal or tumor coverage for allele $allele from $resultdir/${allele}.coverage.txt — check calculate_coverage.sh output format."
    fi

    printf "${allele}\t${normal}\t${tumor}\n" >> $resultdir/coverage.txt
done

log_info "####################################################"
log_info "################### ANNOTATION #####################"
log_info "####################################################"

log_info "Compiling mutations ..."
printf "CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tNormal\tTumor\n" > $resultdir/somatic_mutations.txt
empty_size_mml=$(stat -c%s $resultdir/somatic_mutations.txt)

for allele in $var
do
    grep -v ^# $resultdir/$allele.SNP.vcf >> $resultdir/somatic_mutations.txt
    grep -v ^# $resultdir/$allele.DNP.vcf | awk 'length($4) == 2 && length($5) == 2' >> $resultdir/somatic_mutations.txt
done

if [ $(wc -l < $resultdir/somatic_mutations.txt) -eq 1 ]; then
        log_warn "No somatic mutations found in $tumorbam"
         printf "CHROM\tPOS\tREF\tALT\tFILTER\tNormal_Ref\tNormal_Mut\tTumor_Ref\tTumor_Mut\tTumor_MAF\tUV\tFeature\tamino_acid_change\tvariant_classification\tCDS_pos\taa_pos\tvariant_type\tArtifacts\tZygosity\ttumor_purity\tzcMAF\n" > $resultdir/hla_somatic_mutations.txt
else
    log_info "Correcting position of indels ..."
    run_command Rscript $hlakit/indel_pos_correction.R -o $resultdir -m $resultdir/somatic_mutations.txt -f $hlakit/resources/hla.fasta

    log_info "Running samtools mpileup on normal and tumor bam files to count ref and alt alleles ..."
    run_command bash $hlakit/run_mpileup.sh -o $resultdir -m $resultdir/somatic_mutations.txt -p $hlakit

    log_info "Deduplicate mate pairs in mpileup output ..."
    run_command Rscript $hlakit/mpileup_mate_dedup.R -n $resultdir/normal_mpileupout.txt -m $resultdir/tumor_mapq0_mpileupout.txt -u $resultdir/tumor_mapqnonzero_mpileupout.txt 

    log_info "Counting Ref and Mut reads ..."
    run_command Rscript $hlakit/mpileup_counts.R -s $resultdir/somatic_mutations.txt -n $resultdir/normal_mpileupout.txt -m $resultdir/tumor_mapq0_mpileupout.txt -u $resultdir/tumor_mapqnonzero_mpileupout.txt 

    log_info "Finding UV mutations ..."
    run_command Rscript $hlakit/uv_mutations.R -m $resultdir/somatic_mutations.mpileup_counts.txt -f $hlakit/resources/hla.fasta

    log_info "Adding features to mutation list ..."
    run_command bash $hlakit/add_features.sh -i $resultdir/somatic_mutations.mpileup_counts.UV.txt -b $hlakit/resources/hla_bedfile.bed

    log_info "Annotating mutations ..."
    run_command Rscript $hlakit/annotate_mutations.R -m $resultdir/somatic_mutations.mpileup_counts.UV.features.txt -g $hlakit/resources/hla.gtf -b $hlakit/resources/hla.bed -f $hlakit/resources/hla.fasta -p $tumorpurity

    log_info "Deduplicating mutations ..."
    run_command bash $hlakit/duplicated_mutations.sh -m $resultdir/somatic_mutations.mpileup_counts.UV.features.annotated.txt -r $resultdir

    log_info "Add zygosity-corrected Mutant Allele Frequency ..."
    awk -v FS='\t' -v OFS='\t' -v purity="$tumorpurity" 'NR==FNR {
        zygosity[$1] = $2
        next
    }
    FNR==1 { print $0 "\tZygosity\ttumor_purity\tzcMAF"; next }
    {
        split($1, a, "_")
        gene = a[1] "_" a[2]
        status = (gene in zygosity) ? zygosity[gene] : "NA"
        zcMAF = (status == "Homozygous") ? $16 * 2 : $16

        if (status == "Homozygous" && zcMAF < 0.4 * purity)
            $24 = ($24 == "") ? "artifact_low_tumor_zcMAF" : $24 ",artifact_low_tumor_zcMAF"

        print $0 "\t" status "\t" purity "\t" zcMAF
    }' "$resultdir/zygosity.txt" $resultdir/somatic_mutations.mpileup_counts.UV.features.annotated.dedup.txt > $resultdir/somatic_mutations.mpileup_counts.UV.features.annotated.dedup.zc.txt

    log_info "Creating final mutation list ..."
    awk '{print $1"\t"$2"\t"$4"\t"$5"\t"$7"\t"$12"\t"$13"\t"$14"\t"$15"\t"$16"\t"$17"\t"$18"\t"$19"\t"$20"\t"$21"\t"$22"\t"$23"\t"$24"\t"$25"\t"$26"\t"$27}' FS="\t" $resultdir/somatic_mutations.mpileup_counts.UV.features.annotated.dedup.zc.txt > $resultdir/hla_somatic_mutations.txt
    log_info "Final mutation list: $resultdir/hla_somatic_mutations.txt"

    awk '{gsub(/TRUE/, "T")}1' $resultdir/hla_somatic_mutations.txt > $resultdir/hla_somatic_mutations.tmp.txt
    mv $resultdir/hla_somatic_mutations.tmp.txt $resultdir/hla_somatic_mutations.txt

    if [[ "$shortread" == "yes" ]] || [ "$shortread" = "yes" ] || [[ $shortread == yes ]]; then
        log_info "Filtering mutations covered exclusively by short (sub-full-length) reads ..."
        bash $hlakit/short_readlength_filter.sh -o $resultdir
    fi

    log_info "Creating filtered final mutation list ..."
    awk -F'\t' 'NR==1 || ($5=="PASS" && $18=="")' FS="\t" $resultdir/hla_somatic_mutations.txt > $resultdir/hla_somatic_mutations_filt.txt
    log_info "Final filtered mutation list: $resultdir/hla_somatic_mutations_filt.txt"
fi


log_info "Organizing results ..."
mkdir -p $resultdir/coverage $resultdir/bamfiles/SNP $resultdir/bamfiles/DNP $resultdir/somatic_mutations/SNP $resultdir/somatic_mutations/DNP $resultdir/totrash
mv $resultdir/*coverage.txt $resultdir/coverage/.
mv $resultdir/allelelist_updated.txt $resultdir/bamfiles/.
mv $resultdir/*_subsetted.${reference}.bam* $resultdir/bamfiles/.
mv $resultdir/hla_*SNP.bam* $resultdir/hla_*SNP.MAPQcorrected.bam* $resultdir/bamfiles/SNP/.
mv $resultdir/hla_*DNP.bam* $resultdir/hla_*DNP.MAPQcorrected.bam* $resultdir/bamfiles/DNP/.
mv $resultdir/hla_somatic_mutations*.txt $resultdir/somatic_mutations/.
mv $resultdir/*SNP.vcf $resultdir/somatic_mutations/SNP/.
mv $resultdir/*DNP.vcf $resultdir/somatic_mutations/DNP/.
mv $resultdir/shortread_filter $resultdir/totrash/.
find $resultdir -maxdepth 1 -type f -exec mv {} $resultdir/totrash/ \;
mv $resultdir/totrash/*log $resultdir/.
rm -r $resultdir/totrash

log_info "####################################################"
log_info "###################### LOH #########################"
log_info "####################################################"
loh=$(echo "$loh" | tr -d '\r')
if [[ "$loh" == "yes" ]]; then
    log_info "Analysing Loss of Heterozygosity ..."
    if [ -z "$tumorpurity" ]; then
        log_error "Tumor purity parameter is empty. Cannot analyze Loss of Heterozygosity without Tumor purity."
    else
        if [ -z "$WESnormalcoverage" ] || [ -z "$WEStumorcoverage" ]; then
            log_error "hlakit needs exome-wide coverage for both normal and tumor sample."
        else
            mkdir -p $resultdir/loh
            run_command Rscript $hlakit/loh.R -p $tumorpurity -n $WESnormalcoverage -t $WEStumorcoverage -c $resultdir/coverage/coverage.txt -r $resultdir/loh
            run_command Rscript $hlakit/plot_coverage.R -c $resultdir/coverage/coverage.txt -n $WESnormalcoverage -t $WEStumorcoverage --aibfile $resultdir/loh/aib.txt --lowcov $resultdir/loh/low_cov.txt -r $resultdir/loh
        fi
    fi
fi