# Prefix for references and genetic maps:

# Now pointing to a local resource of the parent project.

# Probably using a public project resource in the future?

refproject="$DX_PROJECT_CONTEXT_ID:/refs"

# available references

refbase_1kg_hg19="/qrefs/1000G/hg19/##.genotypes.qref"

refbase_1kg_hg38lift="/qrefs/1000G/hg38.lifted/chr##.20130502.genotypes.hg38lifted.qref"

refbase_ukb_hg19="/qrefs/UKB/hg19.w20250818/##.ukb23374_v1_cfilter_restored_mac4.phased.w139525_20250818.hg19lifted.qref"

refbase_ukb_hg38="/qrefs/UKB/hg38.w20250818/chr##.ukb23374_v1_cfilter_restored_mac4.phased.w139525_20250818.qref"

refbase_ukb_maf_hg19="/qrefs/UKB/hg19.maf0.0001.w20250818/##.ukb23374_v1_cfilter_restored_mac4.phased.maf0.0001.w139525_20250818.hg19lifted.qref"

refbase_ukb_maf_hg38="/qrefs/UKB/hg38.maf0.0001.w20250818/chr##.ukb23374_v1_cfilter_restored_mac4.phased.maf0.0001.w139525_20250818.qref"

# genetic maps

gmbase_hg19="/1000G/genetic_maps/hg19/genetic_map_hg19_chr##.txt"

gmbase_hg38="/1000G/genetic_maps/hg38/genetic_map_hg38_chr##.txt"

# sample file:

# - either a simple one-column text file with IDs corresponding to the samples in the reference panels

# - or a PLINK .fam file (no header, sample IDs still must be in the first (FID) column (mostly FID and IID are the same))

# - a BGEN .sample file (two header lines are ignored, sample IDs are taken from the second column)

# The sample file is used to mask withdrawn participants from the panel. For this purpose, withdrawals

# have either a sample ID starting with "W" followed by digits or a negative sample ID, i.e. starting

# with "-" followed by digits.

samplefile="${refproject}/qrefs/UKB/samples_w139525"

echo "target: '$target'"

echo "build: '$build'"

echo "reference: '$reference'"

echo "skipPhasing: '$skipPhasing'"

echo "skipImputation: '$skipImputation'"

echo "imputeInfo: '$imputeInfo'"

echo "imputeR2filter: $imputeR2filter"

echo "imputeMAFfilter: $imputeMAFfilter"

echo "outputPhasedFile: '$outputPhasedFile'"

echo "outputUnphased: '$outputUnphased'"

echo "allowRefAltSwap: '$allowRefAltSwap'"

echo "allowStrandFlip: '$allowStrandFlip'"

echo "K: $K"

echo -n "maxChunkMemory: "

if (( $maxChunkMemory <= 0 )); then

maxChunkMemory=0

echo "auto"

else

echo "$maxChunkMemory GiB"

fi

# Download target file into VM.

# To recover the original filename, the output of "dx describe

# "$target" --name" is used.

dx download "$target"

tgt=$(dx describe "$target" --name)

# To report any recognized errors in the correct format in

# $HOME/job_error.json the dx-jobutil-report-error utility

# is used as follows:

#

# dx-jobutil-report-error "My error message"

#

# However, this entire bash script is executed with -e

# when running in the cloud, so any line which returns a nonzero

# exit code will prematurely exit the script; if no error was

# reported in the job_error.json file, then the failure reason

# will be AppInternalError with a generic error message.

# get chromosome number from input file name (it's a relative file name without path already)

chrprefixed=${tgt%%.*} # strip all from the first dot

chr=${chrprefixed/chr/} # remove "chr" -> now, only the number should be left

# determine input format -> will also be the output format

if [[ $tgt == *.bcf ]]; then

fmt=b

ext=".bcf"

else # we only have these two options

fmt=z

ext=".vcf.gz"

fi

# parse imputation options:

unset eagleimpopts

if [[ $skipPhasing == true && $skipImputation == true ]]; then

dx-jobutil-report-error "ERROR: skipImputation and skipPhasing cannot be set together!"

exit 1

fi

if [[ $skipPhasing == true ]]; then

eagleimpopts="$eagleimpopts --skipPhasing"

fi

if [[ $skipImputation == true ]]; then

eagleimpopts="$eagleimpopts --skipImputation"

# override outputPhasedFile option as the phasing output will be required when imputation is skipped

outputPhasedFile=true

fi

unset r2filter

unset maffilter

if [[ $skipImputation != true ]]; then # do imputation

eagleimpopts="$eagleimpopts --imputeInfo $imputeInfo"

if (( $(echo "$imputeR2filter > 0" | bc -l) )); then

eagleimpopts="$eagleimpopts --imputeR2filter $imputeR2filter"

r2filter=1

fi

if (( $(echo "$imputeMAFfilter > 0" | bc -l) )); then

eagleimpopts="$eagleimpopts --imputeMAFfilter $imputeMAFfilter"

maffilter=1

fi

fi

if [[ $skipPhasing != true ]]; then # do phasing

eagleimpopts="$eagleimpopts --K $K"

if [[ $outputPhasedFile == true ]]; then

eagleimpopts="$eagleimpopts --outputPhasedFile"

fi

if [[ $outputPhasedFile == true && $outputUnphased == true ]]; then # available only in combination

eagleimpopts="$eagleimpopts --outputUnphased"

fi

fi

if [[ $allowRefAltSwap == true ]]; then

eagleimpopts="$eagleimpopts --allowRefAltSwap"

fi

if [[ $allowStrandFlip == true ]]; then

eagleimpopts="$eagleimpopts --allowStrandFlip"

fi

# names of final result files

resultadd=."$reference"

if [[ -n "$r2filter" ]]; then

resultadd="$resultadd"."${imputeR2filter}r2"

fi

if [[ -n "$maffilter" ]]; then

resultadd="$resultadd"."${imputeMAFfilter}maf"

fi

tgtbase="${tgt%$ext}"

imputed="$tgtbase$resultadd".imputed$ext

phased="$tgtbase$resultadd".phased$ext

confidences="$tgtbase$resultadd".phased.confidences

varinfo="$tgtbase$resultadd".varinfo

logfile="$tgtbase$resultadd".log

# Everything from here will be logged in a separate log file that

# will be made available as another result file.

{

echo

echo "Starting EagleImp imputation app:"

echo

echo -n "Check if input filename complies with the specification... "

# check if the file is named correctly

check=0

if [[ $chr =~ ^[0-9]+$ ]]; then

if (( $chr > 23 )); then

check=1

fi

elif [[ $chr != "X" && $chr != "X_PAR1" && $chr != "X_nonPAR" && $chr != "X_PAR2" && $chr != "23_PAR1" && $chr != "23_nonPAR" && $chr != "23_PAR2" ]]; then

check=1

fi

if [[ $check != 0 ]]; then

echo "no"

dx-jobutil-report-error "ERROR: Ensure that your target file is named in compliance with the specification, i.e. the chromosome number (1-23) at the beginning separated with a '.'-character. A preceding 'chr' is allowed as well as 'X' instead of the number. 'X' or '23' may also be suffixed with '_PAR1', '_nonPAR' or '_PAR2'."

exit 1

else

echo "yes"

fi

# VM specifications

mem="$(grep MemFree /proc/meminfo | awk '{print $2}')"

disk="$(df -P . | tail -1 | awk '{print $4}')"

echo

echo "VM specifications:"

echo " Available RAM: $mem kiB == $( bc <<< 'scale=3; '$mem' / 1024 / 1024' ) GiB"

echo " Available disk space: $disk kiB == $( bc <<< 'scale=3; '$disk' / 1024 / 1024' ) GiB"

echo -n "Check if the VM is capable to handle the imputation... "

if [[ "$reference" != "1000G" ]] && (( $mem < 120000000 || $disk < 200000000 )); then

echo "no"

dx-jobutil-report-error "ERROR: Ensure that the VM has at least 120 GB of RAM and 200 GB of disk space for UKB imputation."

exit 1

fi

echo "yes"

# take 85% of the available RAM for imputation chunks (and convert to GiB) (auto setting)

mem="$(( $mem * 85 / 100 / 1024 / 1024 ))"

if (( $mem < $maxChunkMemory )); then

echo "WARNING: Your individual setting for maxChunkMemory may exceed the VM's capabilities which could likely result in a system crash!"

fi

# take over individual setting if applicable

if (( $maxChunkMemory > 0 )); then

mem=$maxChunkMemory

fi

# check if filter values are valid

echo

echo -n "Check if filter values are valid... "

if (( $(echo "$imputeR2filter < 0" | bc -l) || $(echo "$imputeR2filter > 1" | bc -l) || $(echo "$imputeMAFfilter < 0" | bc -l) || $(echo "$imputeMAFfilter > 0.5" | bc -l) )); then

echo "no"

dx-jobutil-report-error "ERROR: Invalid filters. Ensure 0 <= R2 <= 1 and 0 <= MAF <= 0.5."

exit 1

else

echo "yes"

fi

# check, if this is chromosome X (without explicit PAR region declaration)

# -> splitting required

splitx=0

if [[ $chr == 23 || $chr == X ]]; then

splitx=1

fi

# choose reference and genetic map according to user selection

if [[ "$build" == "hg19" ]]; then

gmbase="$gmbase_hg19"

if [[ "$reference" == "UKBcomplete" ]]; then

refbase="$refbase_ukb_hg19"

elif [[ "$reference" == "UKBmaf" ]]; then

refbase="$refbase_ukb_maf_hg19"

else

refbase="$refbase_1kg_hg19"

fi

else # hg38

gmbase="$gmbase_hg38"

if [[ "$reference" == "UKBcomplete" ]]; then

refbase="$refbase_ukb_hg38"

elif [[ "$reference" == "UKBmaf" ]]; then

refbase="$refbase_ukb_maf_hg38"

else

refbase="$refbase_1kg_hg38lift"

fi

fi

# choose genetic map file

# special treatment for chrX as we don't seperate the PAR regions for the genetic map

if [[ $splitx == 0 && ( $chr == X* || $chr == 23* ) ]]; then

# a PAR/nonPAR region is directly addressed

if [[ $chr == X* ]]; then

genmap="$refproject${gmbase/\#\#/X}"

else

genmap="$refproject${gmbase/\#\#/23}"

fi

else

genmap="$refproject${gmbase/\#\#/$chr}"

fi

echo

echo -n "Check if genetic map $genmap is available... "

if ! {

dx describe "$genmap" >/dev/null 2>&1

}; then

echo "no"

dx-jobutil-report-error "ERROR: Ensure that the genetic map file is available at the correct location in your project: $genmap"

exit 1

else

echo "yes"

fi

echo -n "Check if sample file $samplefile is available... "

if ! {

dx describe "$samplefile" >/dev/null 2>&1

}; then

echo "no"

dx-jobutil-report-error "ERROR: Ensure that the sample file is available at the correct location in your project: $samplefile"

exit 1

else

echo "yes"

fi

echo -n "Downloading genetic map... "

if ! {

dx download --no-progress $genmap

}; then

dx-jobutil-report-error "ERROR: Failed to download genetic map $genmap"

exit 1

fi

echo "done"

echo -n "Downloading sample file... "

if ! {

dx download --no-progress $samplefile

}; then

dx-jobutil-report-error "ERROR: Failed to download sample file $samplefile"

exit 1

fi

echo "done"

# strip file path from downloaded files as they are located directly in the home folder now

genmap="${genmap##*/}"

samplefile="${samplefile##*/}"

# start timestamp for file processing

procstart=$(date +%s)

# split X, if required

if [[ $splitx == 0 ]]; then

# usual processing of a single file

# reference file (replaces ## by the chromosome number)

ref="$refproject${refbase/\#\#/$chr}"

# process file

process_file "$tgt" "$ref"

# rename output files

imp="${tgt%$ext}".imputed$ext

if [[ -e $imp ]]; then

mv "$imp" "$imputed"

fi

phs="${tgt%$ext}".phased$ext

if [[ -e $phs ]]; then

mv "$phs" "$phased"

fi

conf="${tgt%$ext}".phased.confidences

if [[ -e $conf ]]; then

mv "$conf" "$confidences"

fi

var="${tgt%$ext}".varinfo

mv "$var" "$varinfo"

else

# split X chromosome in PAR/nonPAR regions and process parts separately with final merging

echo -n "Chromosome X requires splitting in PAR/nonPAR regions. Splitting... "

if [[ "$build" == "hg19" ]]; then

regpar1="1-2699520"

regnonpar="2699521-154931043"

regpar2="154931044-"

else

regpar1="1-2781479"

regnonpar="2781480-155701382"

regpar2="155701383-"

fi

if [[ $fmt == b ]]; then

tgt_par1="${chrprefixed}_PAR1.bcf"

tgt_nonpar="${chrprefixed}_nonPAR.bcf"

tgt_par2="${chrprefixed}_PAR2.bcf"

else

tgt_par1="${chrprefixed}_PAR1.vcf.gz"

tgt_nonpar="${chrprefixed}_nonPAR.vcf.gz"

tgt_par2="${chrprefixed}_PAR2.vcf.gz"

fi

if ! {

bcftools index $tgt

bcftools view $tgt -r "23:$regpar1,chr23:$regpar1,X:$regpar1,chrX:$regpar1" -O$fmt -o "$tgt_par1"

bcftools view $tgt -r "23:$regnonpar,chr23:$regnonpar,X:$regnonpar,chrX:$regnonpar" -O$fmt -o "$tgt_nonpar"

bcftools view $tgt -r "23:$regpar2,chr23:$regpar2,X:$regpar2,chrX:$regpar2" -O$fmt -o "$tgt_par2"

}; then

dx-jobutil-report-error "ERROR: Splitting in PAR/nonPAR regions failed."

exit 1

fi

echo "done"

# Processing PAR/nonPAR regions

echo

echo "Processing PAR1..."

ref_par1="$refproject${refbase/\#\#/${chr}_PAR1}"

process_file "$tgt_par1" "$ref_par1"

echo

echo "Processing nonPAR..."

ref_nonpar="$refproject${refbase/\#\#/${chr}_nonPAR}"

process_file "$tgt_nonpar" "$ref_nonpar"

echo

echo "Processing PAR2..."

ref_par2="$refproject${refbase/\#\#/${chr}_PAR2}"

process_file "$tgt_par2" "$ref_par2"

# merging result files

echo

echo "Merging/concatenating result files..."

# imputed

imp_par1="${tgt_par1/PAR1/PAR1.imputed}"

imp_nonpar="${tgt_nonpar/nonPAR/nonPAR.imputed}"

imp_par2="${tgt_par2/PAR2/PAR2.imputed}"

imp_files=""

if [[ -e $imp_par1 ]]; then

imp_files="$imp_files $imp_par1"

fi

if [[ -e $imp_nonpar ]]; then

imp_files="$imp_files $imp_nonpar"

fi

if [[ -e $imp_par2 ]]; then

imp_files="$imp_files $imp_par2"

fi

# phased

phs_par1="${tgt_par1/PAR1/PAR1.phased}"

phs_nonpar="${tgt_nonpar/nonPAR/nonPAR.phased}"

phs_par2="${tgt_par2/PAR2/PAR2.phased}"

phs_files=""

if [[ -e $phs_par1 ]]; then

phs_files="$phs_files $phs_par1"

fi

if [[ -e $phs_nonpar ]]; then

phs_files="$phs_files $phs_nonpar"

fi

if [[ -e $phs_par2 ]]; then

phs_files="$phs_files $phs_par2"

fi

if ! {

if [[ -n "$imp_files" ]]; then

bcftools concat $imp_files -n -o "$imputed"

fi

if [[ -n "$phs_files" ]]; then

bcftools concat $phs_files -n -o "$phased"

fi

}; then

dx-jobutil-report-error "ERROR: Concatenating PAR/nonPAR regions failed."

exit 1

fi

# confidences

shopt -s nullglob # don't throw an error if files referenced with '*' don't exist

for conf in *PAR1.phased.confidences *nonPAR.phased.confidences *PAR2.phased.confidences; do

echo "$conf:" >> "$confidences"

cat "$conf" >> "$confidences"

echo >> "$confidences"

done

# varinfo

header=0

for vinfo in *PAR1.varinfo *nonPAR.varinfo *PAR2.varinfo; do

if [[ header == 0 ]]; then

cat "$vinfo" > "$varinfo"

header=1

else

tail -n +2 "$vinfo" >> "$varinfo"

fi

done

fi # end if splitx...

# end timestamp for file processing

procend=$(date +%s)

proctime=$(( procend - procstart ))

# start timestamp for upload

uplstart=$procend

echo

echo -n "Uploading result files..."

# Upload result files from the VM to the parent project.

# The utility dx-jobutil-add-output is used to format and

# add output variables to the job's output as appropriate for the output

# class.

if [[ -e $imputed ]]; then

imputed_file=$(dx upload "$imputed" --brief)

dx-jobutil-add-output imputed_file "$imputed_file" --class=file

fi

if [[ -e $phased ]]; then

phased_file=$(dx upload "$phased" --brief)

dx-jobutil-add-output phased_file "$phased_file" --class=file

fi

if [[ -e $confidences ]]; then

phasing_confidences=$(dx upload "$confidences" --brief)

dx-jobutil-add-output phasing_confidences "$phasing_confidences" --class=file

fi

varinfo=$(dx upload "$varinfo" --brief)

dx-jobutil-add-output varinfo "$varinfo" --class=file

# end timestamp for upload

uplend=$(date +%s)

upltime=$(( uplend - uplstart ))

echo "done"

echo "Upload time: $upltime seconds"

echo

echo "Total time: $(( proctime + upltime )) seconds"

# log ends here

} 2>&1 | tee -a log

# Whenever an error was thrown, the process will have exited with an error code.

# Unfortunately, the return value here belongs to "tee", but we can access PIPESTATUS.

ret=${PIPESTATUS[0]}

if [[ $ret != 0 ]]; then

# simply exit as the error reason was already set

exit $ret

fi

# This is to make the logfile name conform to the other output files.

mv log "$logfile"

# upload log

log=$(dx upload "$logfile" --brief)

dx-jobutil-add-output log "$log" --class=file

exit 0

# We got tgt and ref as parameters, so I make them local here to shadow the global context.

local tgt="$1"

local ref="$2"

echo

echo -n "Check if reference file $ref is available... "

if ! {

dx describe "$ref" >/dev/null 2>&1

}; then

echo "no"

dx-jobutil-report-error "ERROR: Ensure that the reference file is available at the correct location in your project: $ref"

exit 1

else

echo "yes"

fi

# start timestamp for reference download

dwnlstart=$(date +%s)

echo -n "Downloading reference... "

if ! {

dx download --no-progress $ref

}; then

dx-jobutil-report-error "ERROR: Failed to download reference $ref"

exit 1

fi

echo "done"

# end timestamp

dwnlend=$(date +%s)

downloadtime=$(( dwnlend - dwnlstart ))

echo "Download time: $downloadtime seconds"

# strip file path from just downloaded file as it is located directly in the home folder now

ref="${ref##*/}"

# start timestamp for index

idxstart=$(date +%s)

echo

echo -n "Run bcftools index $tgt... "

if ! {

bcftools index $tgt

}; then

echo "Index failed."

dx-jobutil-report-error "ERROR: bcftools index $tgt failed."

exit 1

fi

echo "done"

# end timestamp for index

idxend=$(date +%s)

idxtime=$(( idxend - idxstart ))

echo "Indexing time: $idxtime seconds"

# EagleImp command

cmd="eagleimp --target $tgt --ref $ref --geneticMap $genmap --excludeWithdrawals $samplefile -O$fmt --maxChunkMemory $mem $eagleimpopts"

# start timestamp for eagleimp

eagleimpstart=$(date +%s)

echo

echo "Run EagleImp..."

echo "$cmd"

if ! {

$cmd 2>&1

}; then

dx-jobutil-report-error "ERROR: Eagleimp failed. Command: $cmd"

exit 1

fi

# end timestamp for eagleimp

eagleimpend=$(date +%s)

eagleimptime=$(( eagleimpend - eagleimpstart ))

echo

echo "EagleImp time: $eagleimptime seconds"

echo "Total processing time for file $tgt: $(( downloadtime + idxtime + eagleimptime )) seconds"