#!/bin/bash
set -e

# requires: bwa, samtools, bcftools
# if using modules:
source /mnt/software/Modules/current/init/bash
module load bwa samtools bcftools

if [ $# -lt 2 ]; then
        echo "Usage: $0 ref /path/to/laa_job"
        echo "Output: Per barcode alignment to reference (*.aln.bam) and variant calls (*.vcf)"
        exit 1
fi

REF="$1"
FQTGZ="$2/tasks/pbcoretools.tasks.split_laa_fastq-0/consensus_fastq.tar.gz"
FQZ="$2/tasks/pbcoretools.tasks.split_laa_fastq-0/consensus_fastq.zip"

# make sure the reference is indexed
samtools faidx "${REF}"
bwa index "${REF}"

# unzip fastqs
if [ -e "${FQTGZ}" ]; then
        tar zxvf "${FQTGZ}"
elif [ -e "${FQZ}" ]; then
        unzip -j "${FQZ}"
else
        echo "LAA job not found at $2"
        exit 1
fi

# align to reference and call variants
for fastq in consensus_fastq.*.fastq
do
        bwa mem "${REF}" "${fastq}" | samtools view -Sb - | samtools sort - > "${fastq%.*}.aln.bam"
        samtools index ${fastq%.*}.aln.bam
        samtools mpileup -uf "${REF}" "${fastq%.*}.aln.bam" | bcftools call -mv --ploidy 1 - > "${fastq%.*}.vcf"
done

module remove bwa samtools bcftools
set +e