14_genotypeDB.sh

#!/bin/bash

#############
## written 2025.04.08 RMPeery
## USE: ./genotypeDB.sh
## Note: this is a loop script to combine gVCFs per interval, 
## with one interval per loop then call genotypes across all samples
#############

set -eou pipefail

# change to match your project
DIR=/home/user/path
output_vcfs="$DIR/vcfs/intervalVCFs_splitBed"
reference="$DIR/refs/Phyci1_AssemblyScaffolds.fasta"
interval_dir="$DIR/refs/beds"
gvcfs="$DIR/gvcf"
interval_dbs="$DIR/vcfs/interval_dbs_splitBed"

# safely make directories jic
mkdir -p "$output_vcfs"
mkdir -p "$interval_dbs"

# name the dataset
DATASET=Pcin

# put gvcfs into an array
gvcf_in=("$gvcfs"/*.g.vcf.gz)

# make the genotype call command
# input from ChatGPT3 on this
cmd=""
for file in "${gvcf_in[@]}"; do
    cmd+="--variant $file "
done

# run intervals sequentially
for bed in "$interval_dir"/interval_batch*.bed; do
    bedfile=$(basename "$bed")
    db_path="$interval_dbs/${bedfile}_db"
    out_vcf="$output_vcfs/${bedfile}.vcf.gz"
    # write out to check inputs are found
    echo "Processing: $bedfile"
    echo "Creating GenomicsDB: $db_path"
    echo "Output VCF: $out_vcf"

    # remove db if it exists - otherwise the run fails
    rm -rf "$db_path"

    # run gatk GenomicsDBImport with extra java options for memory
    gatk --java-options "-Xmx20g -Xms12g" GenomicsDBImport \
        $cmd \
        --genomicsdb-workspace-path "$db_path" \
        --reader-threads 2 \
        -L "$bed"

    # run gatk GenotypeGVCFs with extra java options for memory
    gatk --java-options "-Xmx32g -DGATK_STACKTRACE_ON_USER_EXCEPTION=true" GenotypeGVCFs \
        -R "$reference" \
        -V "gendb://$db_path" \
        -O "$out_vcf"

    echo "Finished processing: $bedfile"
done

echo "All interval genotyping completed."