#!/bin/bash
#PBS -q queue_name
#PBS -l walltime=240:00:00
#PBS -l nodes=1,flags=allprocs
#PBS -j oe
#PBS -N arrow_polish

bin=/path/to/local/opt/miniconda2/bin
path=/path/to/working/directory
dest=$path/output_directory
ref=$path/assembly.ctg.fa
name=asm_ctg
threads=24

samples=(
    $path/pacbio_sequel_ii_clr.run1.subreads.bam
    $path/pacbio_sequel_ii_clr.run2.subreads.bam
)

cd $dest

#
# Align and sort the raw reads using pbmm2, which is PacBio's wrapper for minimap2
#
for sample in "${samples[@]}"
do
    out=`basename $sample | sed -E 's/.subreads.bam//'`
    $bin/pbmm2 align -j $threads --sort $ref $sample > ${dest}/${out}.bam
    mergelist="${dest}/${out}.bam "$mergelist
done

#
# Merge the aligned files.
#
samtools merge --threads 8 $dest/${name}.bam $mergelist

#
# Remove the unmerged files.
#
for sample in $mergelist
do
    rm $sample
done

#
# Index the PacBio-specific data in the BAM file
#
$bin/pbindex $dest/${name}.bam
samtools faidx $ref

#
# Run Arrow variant Caller
#
$bin/variantCaller \
    --algorithm arrow \
    --numWorkers $threads \
    --referenceFilename $ref \
    -o $dest/${name}.arrow.fa \
    -o $dest/${name}.arrow.fq \
    -o $dest/${name}.arrow.gff \
    $dest/${name}.bam
