## --------------------------------------------------------------
##  Figure with 10 bar-plots (2 × 5) – one page only
##  Produced by Sven BUERKI, October 2025
## --------------------------------------------------------------

## ---- 1. Load & filter data ------------------------------------
phytop <- read.csv("ASTRAL_renamed_final.tre.info.tsv", sep = "\t")

## keep only nodes that have a label (ILS/IH nodes)
phytop_f <- subset(phytop, Node_Label_ASTRAL_Tree != "")

## ---- 2. Keep at most 10 panels --------------------------------
n_plots <- min(10, nrow(phytop_f))
phytop_f <- phytop_f[1:n_plots, ]

## ---- 3. Pre-compute legend text (one block per panel) ----------
legends <- with(phytop_f, paste0(
  "n=", ceiling(n), "\n",
  "P-value=", format(p_value, scientific = TRUE), "\n",
  "ILS-e=", round(100 * ILS_explain, 1), "%\n",
  "IH-e=",  round(100 * IH_explain,  1), "%\n",
  "ILS-i=", round(100 * ILS_index,   1), "%\n",
  "IH-i=",  round(100 * IH_index,    1), "%"
))

## ---- 4. Open a *large* graphics device (optional but safe) -----
# Choose ONE of the three lines:
pdf("Phytop_quartet_10panels_nodeswithILS&IH.pdf", width = 14, height = 8)   # vector PDF

## ---- 5. Layout & margins ---------------------------------------
par(
  mfrow = c(2, 5),          # 2 rows × 5 columns = 10 panels
  mar   = c(4.5, 4.5, 3, 7),# bottom, left, top, right (7 lines for legend)
  oma   = c(0, 0, 3, 0),    # outer margin for optional overall title
  xpd   = TRUE              # allow drawing outside plot region
)

## ---- 6. Plot loop -----------------------------------------------
panel_letters <- LETTERS[1:n_plots]   # A, B, … J

for (i in seq_len(n_plots)) {
  
  ## quartet scores
  q <- c(phytop_f$q1[i], phytop_f$q2[i], phytop_f$q3[i])
  
  ## ---- barplot ------------------------------------------------
  barplot(q,
          col       = c("blue", "orange", "grey"),
          names.arg = c("q1", "q2", "q3"),
          ylab      = "Quartet score",
          ylim      = c(0, 1),
          las       = 1)                 # vertical axis labels
  
  ## ---- panel label (A–J) --------------------------------------
  mtext(panel_letters[i],
        side = 3, adj = 0, line = 0.5,
        cex = 1.4, font = 2)           # bold
  
  ## ---- legend (outside, right side) ---------------------------
  legend_text <- unlist(strsplit(legends[i], "\n"))
  
  legend("topright",
         legend = legend_text,
         bty    = "n",          # no box
         inset  = c(-0.38, 0), # push ~38 % of plot width to the right
         cex    = 0.95,
         text.col = "black")
}

## ---- 7. Close device -------------
dev.off() 




