## ----knitr_options, include = FALSE-------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    warning = FALSE,
    message = FALSE
)

## ----libraries----------------------------------------------------------------
library(ape)
library(dplyr)
library(phylobar)
library(readr)
library(tibble)
library(tidyr)
library(seriation)
library(zen4R)

## ----read_in_data-------------------------------------------------------------
download_zenodo("10.5281/zenodo.17477876", tempdir())
all_data <- readRDS(file.path(tempdir(), "/su-2020.rds"))
cell_counts <- all_data$cell_counts
metadata <- all_data$metadata

## ----compositions-------------------------------------------------------------
x <- cell_counts |>
    dplyr::count(patient_id, majority_voting) |>
    pivot_wider(
        names_from = majority_voting,
        values_from = n,
        values_fill = 0
    ) |>
    column_to_rownames("patient_id")
x <- x / rowSums(x)

## ----make_phylo_from_igraph---------------------------------------------------
edges_text <- "source,target
Hematopoietic stem cell,Common myeloid progenitor
Hematopoietic stem cell,Common lymphoid progenitor
Common myeloid progenitor,Megakaryocyte
Megakaryocyte,Platelet
Common myeloid progenitor,RBC
Common myeloid progenitor,Myeloblast
Myeloblast,Monocyte
Monocyte,CD14 Monocyte
Monocyte,CD16 Monocyte
Myeloblast,SC & Eosinophil
Myeloblast,DC
Myeloblast,pDC
Common lymphoid progenitor,NK
Common lymphoid progenitor,Small lymphocyte
Small lymphocyte,T lymphocyte
T lymphocyte,CD4 T lymphocyte
T lymphocyte,CD8 T lymphocyte
CD4 T lymphocyte,CD4 T
CD4 T lymphocyte,CD4m T
CD4 T lymphocyte,CD4n T
CD8 T lymphocyte,CD8m T
CD8 T lymphocyte,CD8eff T
T lymphocyte,gd T
Small lymphocyte,B lymphocyte
B lymphocyte,B
B lymphocyte,IgA PB
B lymphocyte,IgG PB"

edge <- read_csv(edges_text)
tips <- setdiff(edge$target, edge$source)

# preparation for phylo construction
internal_nodes <- setdiff(unlist(edge), tips)
node_order <- c(tips, internal_nodes)
ix <- setNames(seq_along(node_order), node_order)

# create phylo
tree <- list(
    edge = matrix(c(ix[edge$source], ix[edge$target]), ncol = 2),
    Nnode = length(internal_nodes),
    tip.label = node_order[seq_along(tips)],
    node.label = node_order[seq(length(tips) + 1, length(ix))]
)
class(tree) <- "phylo"

## ----patient-order------------------------------------------------------------
md <- metadata |>
    distinct(clinical_id, severity) |>
    mutate(severity = factor(
        severity, levels = c("mild", "moderate", "severe")
    )) |>
    arrange(severity)

patient_order <- md %>%
    split(.$severity) |>
    lapply(\(df) {
        ids <- df$clinical_id
        m <- x[ids, , drop = FALSE]
        ids[hclust(dist(m))$order]
    }) |>
    unlist(use.names = FALSE)

## ----phylobar-plot------------------------------------------------------------
phylobar(x[patient_order, ], tree)

## ----session------------------------------------------------------------------
sessionInfo()