#!/usr/bin/env bash
# Rebuilds the example tree in the phylogenetic tree guide.
# Nine complete primate mitochondrial genomes from NCBI RefSeq (accessions.tsv),
# aligned with MAFFT 7.526 and analysed with IQ-TREE 3.1.1: ModelFinder,
# 1000 ultrafast bootstrap and 1000 SH-aLRT replicates, rooted on rhesus macaque.
#
# Usage: MAFFT=/path/to/mafft IQTREE=/path/to/iqtree3 bash build_tree.sh
set -euo pipefail
cd "$(dirname "$0")"
MAFFT="${MAFFT:-mafft}"
IQTREE="${IQTREE:-iqtree3}"

if [ ! -s primate_mtdna.fasta ]; then
  ids=$(tail -n +2 accessions.tsv | cut -f2 | paste -sd, -)
  curl -fsS "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=${ids}&rettype=fasta&retmode=text" -o raw.fasta
  # Replace each RefSeq header with the short label from accessions.tsv.
  awk -F'\t' 'NR==FNR { if (FNR > 1) label[$2] = $1; next }
    /^>/ { split(substr($0, 2), f, " "); print ">" label[f[1]]; next }
    { print }' accessions.tsv raw.fasta > primate_mtdna.fasta
  rm raw.fasta
fi

"$MAFFT" --auto --thread 4 primate_mtdna.fasta > primate_mtdna.aln.fasta

rm -rf iqtree && mkdir iqtree
"$IQTREE" -s primate_mtdna.aln.fasta -m MFP -B 1000 -alrt 1000 \
  -o Rhesus_macaque --seed 12345 -T 4 --prefix iqtree/primate_mtdna