Phylogenetics & Tree Building

# MAFFT — fastest, most used
mafft --auto sequences.fasta > aligned.fasta

# MAFFT — high accuracy for <200 sequences
mafft --maxiterate 1000 --localpair sequences.fasta > aligned.fasta

# MUSCLE v5 — newer, competitive accuracy
muscle -align sequences.fasta -output aligned.fasta

# PRANK — for codon-aware alignment (best for Ks/dN/dS)
prank -d=sequences.fasta -o=aligned -codon -F

# trimAl — automated trimming (recommended)
trimal -in aligned.fasta -out trimmed.fasta -automated1

# trimAl — gap threshold (remove columns with >50% gaps)
trimal -in aligned.fasta -out trimmed.fasta -gt 0.5

# Gblocks — stricter trimming
Gblocks aligned.fasta -t=p -b4=5 -b5=h -e=-gb

# ClipKIT — newer, smart trimming
clipkit aligned.fasta -o trimmed.fasta -m kpic-smart-gap

# IQ-TREE's ModelFinder (built-in, best option)
iqtree2 -s trimmed.fasta -m MFP   # MFP = ModelFinder Plus

# Standalone ModelTest-NG
modeltest-ng -i trimmed.fasta -t ml -p 8

# Output example:
# Best model: GTR+F+I+G4 (DNA) or LG+F+G4 (protein)
# BIC score: 45678.234

# Standard analysis: model selection + tree + bootstrap
iqtree2 -s trimmed.fasta \
  -m MFP \
  -bb 1000 \
  -alrt 1000 \
  -nt AUTO \
  --prefix my_tree

# Partition analysis (different genes, different models)
# Create partition file: partition.nex
iqtree2 -s concatenated.fasta \
  -p partition.nex \
  -m MFP \
  -bb 1000 \
  -nt AUTO

# Constrained tree search (test topology)
iqtree2 -s trimmed.fasta -m MFP -g constraint_tree.nwk -bb 1000

# Key output files:
# my_tree.treefile — Newick tree (open in FigTree/iTOL)
# my_tree.iqtree   — full log with model, scores, support
# my_tree.contree  — consensus tree with bootstrap values

# Standard ML tree with bootstrap
raxml-ng --all \
  --msa trimmed.fasta \
  --model GTR+G \
  --bs-trees 200 \
  --threads auto{16} \
  --prefix raxml_tree

# With auto model selection (use modeltest-ng first)
raxml-ng --all --msa trimmed.fasta --model LG+G8+F --bs-trees 200

# Check alignment (sanity check before running)
raxml-ng --check --msa trimmed.fasta --model GTR+G

# DNA (GTR model)
FastTree -gtr -nt trimmed.fasta > tree.nwk

# Protein (default JTT+CAT model)
FastTree trimmed_protein.fasta > tree.nwk

# With local support values (more reliable than default SH)
FastTree -gtr -nt -gamma trimmed.fasta > tree.nwk

# Create a NEXUS file or add MrBayes block to your alignment
begin mrbayes;
  set autoclose=yes nowarn=yes;
  lset nst=6 rates=invgamma;    [GTR+I+G]
  prset brlenspr=unconstrained:exp(10.0);
  mcmc ngen=5000000 samplefreq=1000
       nchains=4 nruns=2
       printfreq=10000;
  sump burnin=1250;
  sumt burnin=1250 contype=allcompat;
end;

# Step 1: Build gene trees for single-copy orthologs
# (From OrthoFinder single-copy orthologs)
for fasta in orthologs/*.fasta; do
  name=$(basename "$fasta" .fasta)
  mafft --auto "$fasta" > "aligned/${name}.aln"
  trimal -in "aligned/${name}.aln" -out "trimmed/${name}.trim" -automated1
  iqtree2 -s "trimmed/${name}.trim" -m MFP -bb 1000 -nt 2 --prefix "trees/${name}"
done

# Step 2: Concatenate all gene trees
cat trees/*.treefile > all_gene_trees.nwk

# Step 3: Run ASTRAL
java -jar astral.jar \
  -i all_gene_trees.nwk \
  -o species_tree.nwk \
  2> astral.log

# Step 4: With branch support (local posterior probability)
java -jar astral.jar -i all_gene_trees.nwk -o species_tree.nwk -t 2

# ggtree — ggplot2-based tree visualization
library(ggtree)
library(treeio)

tree <- read.tree("species_tree.nwk")

# Basic tree
ggtree(tree) +
  geom_tiplab(size=3) +
  geom_nodelab(aes(label=label), size=2, hjust=-0.1) +
  theme_tree2()

# Circular tree with bootstrap colors
ggtree(tree, layout="circular") +
  geom_tiplab2(size=2.5) +
  geom_nodepoint(aes(color=as.numeric(label) > 95), size=2) +
  scale_color_manual(values=c("grey70","red"), labels=c("<95%","≥95%")) +
  theme(legend.position="bottom")

# With heatmap annotation
library(ggtreeExtra)
p <- ggtree(tree) + geom_tiplab()
gheatmap(p, gene_count_matrix, colnames_angle=90, font.size=2)