#!/usr/bin/env bash
# Re-download the sihek genome assembly that is excluded from git.
#
# GCA_033439825.1 (ASM3343982v1) -- Todiramphus cinnamominus, Iridian Genomes,
# released 2023-11-08. The only Guam kingfisher assembly in existence.
# Illumina/SPAdes, scaffold level: 1.16 Gb, 138,582 scaffolds, N50 60,279 bp,
# longest scaffold 590 kb, 0 annotated genes.
#
# Downloads ~353 MB, verifies the MD5 against NCBI's own manifest, then
# decompresses to ~1.1 GB and builds a .fai index.

set -euo pipefail

ACC="GCA_033439825.1_ASM3343982v1"
BASE="https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/033/439/825/${ACC}"
DEST="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)/genome"

mkdir -p "$DEST"
cd "$DEST"

echo "==> metadata"
curl -sSf -O "${BASE}/md5checksums.txt" \
         -O "${BASE}/${ACC}_assembly_report.txt" \
         -O "${BASE}/${ACC}_assembly_stats.txt"

echo "==> genomic FASTA (~353 MB)"
curl -f --progress-bar -O "${BASE}/${ACC}_genomic.fna.gz"

echo "==> verifying MD5 against NCBI's manifest"
expected=$(awk -v f="./${ACC}_genomic.fna.gz" '$2 == f {print $1}' md5checksums.txt)
if [ -z "$expected" ]; then
  echo "FAIL: no checksum for ${ACC}_genomic.fna.gz in md5checksums.txt" >&2
  exit 1
fi
if command -v md5 >/dev/null 2>&1; then
  actual=$(md5 -q "${ACC}_genomic.fna.gz")          # macOS
else
  actual=$(md5sum "${ACC}_genomic.fna.gz" | cut -d' ' -f1)   # GNU
fi
if [ "$expected" != "$actual" ]; then
  echo "FAIL: MD5 mismatch (expected $expected, got $actual)" >&2
  exit 1
fi
echo "    ok  $actual"

echo "==> decompressing (~1.1 GB)"
gunzip -kf "${ACC}_genomic.fna.gz"

echo "==> indexing"
if command -v samtools >/dev/null 2>&1; then
  samtools faidx "${ACC}_genomic.fna"
else
  python3 -c "
from pyfaidx import Faidx
Faidx('${ACC}_genomic.fna')
" 2>/dev/null || {
    echo "    note: install samtools or 'pip3 install --user pyfaidx' to build the .fai index"
    exit 0
  }
fi

echo "==> done"
ls -lh "$DEST"