# Stop on errors.
set -uex

# Reference accession numbers.
ACC=AF086833

# Create the directory for reference file.
mkdir -p refs

# The name of the reference.
REF=refs/${ACC}.fa

# The name of the BAM file
BAM=align.bam

# Obtain the reference genome.
bio fetch -format fasta $ACC > $REF

# Create a bwa index for the reference.
bwa index $REF 1>>log.txt 2>> log.txt

# Create a samtools index for the reference used by IGV.
samtools faidx $REF

# This is the read naming generated by dwgsim.
R1=read1.fq
R2=read2.fq

# Simulate reads from the reference file.
pywgsim -N 10000 $REF $R1 $R2 > simulated.gff

# Generate the alignment from the data simulated above.
bwa mem $REF $R1 $R2 2>> log.txt | samtools sort > $BAM 2>> log.txt

# Index the alignment
samtools index $BAM

# Freebayes with at least 10 reads supporting a variant.
freebayes -C 10 -f $REF $BAM > variants.vcf