# Stop on errors. Show all commands.
set -eux 

# Set the reference genome accession number.
ACC=AF086833

# Set the SRR run number.
SRR=SRR1972739

# Directory for various reference information.
mkdir -p refs

# The reference genome in GENBANK format.
GB=refs/$ACC.gb

# The reference genome in FASTA format.
REF=refs/$ACC.fa

# The reference genome in GFF (Gene Feature Format).
GFF=refs/$ACC.gff

# Get the reference sequence in GenBank format.
efetch -db=nuccore -format=gb -id=$ACC > $GB

# Reformat GenBank to FASTA.
cat $GB | readseq -p -f fasta > $REF 

# Reformat GenBank to GFF. Keep only lines that match CDS.
cat $GB | readseq -p -f gff | grep CDS > $GFF 

# Index reference for the bwa aligner.
bwa index $REF

# Index the reference genome so that it can be loaded into IGV.
samtools faidx $REF

# The directory that stores the sequencing reads
mkdir -p reads

# Get sequencing data from a SRR number.
fastq-dump -X 10000 --split-files -O reads $SRR

# Files names that store read pairs
R1=reads/${SRR}_1.fastq
R2=reads/${SRR}_2.fastq

# The alignment file name.
BAM=$SRR.bam

# All genotype likelihoods.
ALL=$SRR.all.vcf

# Only the variants from the genotype likelhoods.
VCF=$SRR.variants.vcf

# Align and generate a sorted BAM file.
bwa mem $REF $R1 $R2 | samtools sort > $BAM

# Index the BAM file.
samtools index $BAM

# Compute all genotypes with bcftools.
bcftools mpileup -Ou -f $REF $BAM > $ALL

# From all genotypes call only the variants.
cat $ALL | bcftools call --ploidy 1 -mv -Ov -o $VCF 

# Usually it is a waste of space to keep the file with all genoytpes
# We can pipe the two into one call like so.
#bcftools mpileup -Ou -f $REF $BAM | bcftools call -mv -Ov -o $VCF