# Show the commands as are executed.
set -uex

# Create a directory to store the reference and all indices.
mkdir -p refs

# The accession number.
ACC=AF086833

# The reference genome we align to.
REF=refs/REFERENCE.fa

# The annotation file for the reference genome.
GFF=refs/ANNOTATIONS.fa

# The genome file mutated from the reference.
GENOME=GENOME.fa

# The read alignment file.
BAM=align.bam

# The semi-global alignment of the genome to the reference.
GLOBAL=global.bam

# Error rate for the simulation.
ERR=0.05

# How many reads to generate.
N=2500

#
# Install bio with:
#
# pip install bio --upgrade
#
# Obtain the data save the first 3kb of it as our reference FASTA.
bio $ACC --fetch --fasta --end 3kb > $REF

# We create the modified genome here. We want to delete two bases from 1980 to 2000.
# We are replacing a sequence pattern with a shorter version of the same (basically deleting 2 bases).
# You are welcome to experiment with your changes.
cat $REF | biosed -filter -target AGGGTGGACAACAGAAGAACA -replace AGGGTGGACAACAGAAGAA > $GENOME

# Read names
R1=read1.fq
R2=read2.fq

# Generate simulated reads with no other mutations.
wgsim -N $N -e $ERR -r 0 $GENOME $R1 $R2

# Index the reference genome with bwa index.
bwa index $REF

# READ tagging means adding more information into each alignment.
TAG="@RG\tID:variants\tSM:${ACC}\tLB:simulation"

# Align the reads and generate athe BAM file.
bwa mem -R $TAG $REF $R1 $R2 | samtools sort > $BAM

# We also align the genome to the reference.
bwa mem $REF $GENOME | samtools sort > $GLOBAL

# Create the indexs for the BAM files.
samtools index $BAM
samtools index $GLOBAL

# Variant call with bcftools.
bcftools mpileup -Ou -f $REF $BAM | bcftools call -mv -Ov -o bcftools.vcf

# Variant call with freebayes. Filter at 0.1 pval.
freebayes -P 0.1 -f $REF $BAM > freebayes.vcf