## ----echo=FALSE, results="hide", message=FALSE--------------------------------
knitr::opts_chunk$set(error = FALSE, message = FALSE, warning = FALSE)
library(BiocStyle)

## ----message=FALSE, warning=FALSE---------------------------------------------
library(dandelionR)
library(scRepertoire)
library(scater)

## ----warning = FALSE----------------------------------------------------------
data(sce_vdj)

## -----------------------------------------------------------------------------
set.seed(100)

## ----warning = FALSE----------------------------------------------------------
sce_vdj <- setupVdjPseudobulk(sce_vdj,
    already.productive = FALSE,
    allowed_chain_status = c(
        "Single pair", "Extra pair",
        "Extra pair-exception", "Orphan VDJ",
        "Orphan VDJ-exception"
    )
)

## ----warning = FALSE----------------------------------------------------------
plotUMAP(sce_vdj, color_by = "anno_lvl_2_final_clean")

## ----warning = FALSE----------------------------------------------------------
library(miloR)
traj_milo <- Milo(sce_vdj)
milo_object <- buildGraph(traj_milo, k = 50, d = 20, reduced.dim = "X_scvi")
milo_object <- makeNhoods(milo_object, reduced_dims = "X_scvi", d = 20)

## ----warning = FALSE----------------------------------------------------------
milo_object <- miloUmap(milo_object)

## -----------------------------------------------------------------------------
plotUMAP(milo_object,
    color_by = "anno_lvl_2_final_clean",
    dimred = "UMAP_knngraph"
)

## -----------------------------------------------------------------------------
pb.milo <- vdjPseudobulk(milo_object, col_to_take = "anno_lvl_2_final_clean")

pb.milo <- runPCA(pb.milo, assay.type = "Feature_space")

## -----------------------------------------------------------------------------
plotPCA(pb.milo, color_by = "anno_lvl_2_final_clean")

## -----------------------------------------------------------------------------
library(SingleCellExperiment)

# extract the PCA matrix
pca <- t(as.matrix(reducedDim(pb.milo, type = "PCA")))
# define the CD8 terminal cell as the top-most cell and CD4 terminal cell as
# the bottom-most cell
branch.tips <- c(which.max(pca[2, ]), which.min(pca[2, ]))
names(branch.tips) <- c("CD8+T", "CD4+T")
# define the start of our trajectory as the right-most cell
root <- which.max(pca[1, ])

## ----warning = FALSE----------------------------------------------------------
library(destiny)
# Run diffusion map on the PCA
feature_space <- t(assay(pb.milo, "Feature_space"))
dm <- DiffusionMap(as.matrix(feature_space), n_pcs = 50, n_eigs = 10)

## -----------------------------------------------------------------------------
dif.pse <- DPT(dm, tips = c(root, branch.tips), w_width = 0.1)

## ----message=FALSE------------------------------------------------------------
# the root is automatically called DPT + index of the root cell
DPTroot <- paste0("DPT", root)
# store pseudotime in milo object
pb.milo$pseudotime <- dif.pse[[DPTroot]]
# set the colours for pseudotime
pal <- colorRampPalette(rev((RColorBrewer::brewer.pal(9, "RdYlBu"))))(255)
plotPCA(pb.milo, color_by = "pseudotime") +
    scale_colour_gradientn(colours = pal)

## -----------------------------------------------------------------------------
pb.milo <- markovProbability(
    milo = pb.milo,
    diffusionmap = dm,
    terminal_state = branch.tips,
    root_cell = root,
    pseudotime_key = "pseudotime"
)

## ----message=FALSE------------------------------------------------------------
plotPCA(pb.milo, color_by = "CD8+T") + scale_color_gradientn(colors = pal)
plotPCA(pb.milo, color_by = "CD4+T") + scale_color_gradientn(colors = pal)

## -----------------------------------------------------------------------------
cdata <- projectPseudotimeToCell(milo_object, pb.milo, value_key = c("pseudotime", "CD8+T", "CD4+T"))

## ----message=FALSE------------------------------------------------------------
plotUMAP(cdata, color_by = "anno_lvl_2_final_clean", dimred = "UMAP_knngraph")
plotUMAP(cdata, color_by = "pseudotime", dimred = "UMAP_knngraph") +
    scale_color_gradientn(colors = pal)
plotUMAP(cdata, color_by = "CD4+T", dimred = "UMAP_knngraph") +
    scale_color_gradientn(colors = pal)
plotUMAP(cdata, color_by = "CD8+T", dimred = "UMAP_knngraph") +
    scale_color_gradientn(colors = pal)

## ----warning = FALSE----------------------------------------------------------
sessionInfo()