# Stop on errors. Show the commands.
set -uex

# NCBI Bioproject ID.
ID=PRJNA257197

# NCBI accession number for the reference genome.
ACC=AF086833

# The number of samples to process.
N=5

# The directory that stores the reference files.
mkdir -p refs

# The reference data in FASTA format.
REF=refs/${ACC}.fa

# Obtain the reference genome.
efetch -db nuccore -id $ACC -format fasta > $REF

# Index the FASTA file for bwa.
bwa index $REF 2>> log.txt

# Index the FASTA file for IGV.
samtools faidx $REF >> log.txt

# Obtain the SRR run number for the data deposited into the BioProject
esearch -db sra -query $ID | efetch -format runinfo > runinfo.csv

# Keep the top N SRR numbers for now.
cat runinfo.csv  | cut -f 1 -d , | grep SRR | head -$N  > srr.txt

# Make a directory for the sequencing data.
mkdir -p reads

# Limit each dataset to these many reads.
LIMIT=100000

# Download the sequencing data in parallel.
cat srr.txt | parallel "fastq-dump -X $LIMIT -O reads --split-files {} >> log.txt"

# Generate alignments
mkdir -p bam

# Run the bwa mem processes in parallel for all samples.
cat srr.txt | parallel "bwa mem $REF reads/{}_1.fastq reads/{}_2.fastq 2>> log.txt | samtools sort > bam/{}.bam 2>> log.txt"

# Index the BAM files.
cat srr.txt | parallel "samtools index bam/{}.bam 2>> log.txt"

# Call variants on all BAM files at the same time.
bcftools mpileup -Ou -f $REF  bam/*.bam 2>> log.txt | bcftools call --ploidy 1 -vm -Ou |  bcftools norm -Ov -f $REF -d all - > variants.vcf