# Stop on errors.
set -uxeo pipefail

# Set the reference genome accession number.
ACC=AF086833
#ACC=AF086833

# Set the SRR run number.
SRR=SRR1972739
#SRR=SRR1553500

# Get the reference genome.
mkdir -p refs

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

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

# Stores the reference genome in Gene Feature Format.
GFF=refs/$ACC.gff

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

# Reformat GenBank as FASTA
cat $GB | readseq -p -f fasta > $REF 2>> log.txt

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

# Index reference for the aligner.
bwa index $REF 2>> log.txt

# Index the reference genome for IGV.
samtools faidx $REF

# The directory that stores the reads
mkdir -p reads

# Get data from an ebola sequencing run.
fastq-dump -X 10000 --split-files -O reads $SRR

# The name of the read pairs
R1=reads/${SRR}_1.fastq
R2=reads/${SRR}_2.fastq

# This will be the resulting alignment file.
BAM=$SRR.bam
VCF=$SRR.vcf

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

# Index the BAM file.
samtools index $BAM

# Call snps with bcftools
bcftools mpileup -Ou -f $REF $BAM 2>> log.txt | bcftools call -mv -Ov -o $VCF 2>> log.txt