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

## ----eval = FALSE-------------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("dandelionR")

## ----eval = FALSE-------------------------------------------------------------
# if (!requireNamespace("devtools", quietly = TRUE)) {
#     install.packages("devtools")
# }
# devtools::install_github("tuonglab/dandelionR", dependencies = TRUE)

## ----eval = FALSE-------------------------------------------------------------
# # only for the tutorial
# if (!requireNamespace("scater", quietly = TRUE)) {
#     BiocManager::install("scater")
# }
# if (!requireNamespace("scRepertoire", quietly = TRUE)) {
#     BiocManager::install("scRepertoire")
# }
# # or
# devtools::install_github("ncborcherding/scRepertoire")

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

## -----------------------------------------------------------------------------
data(demo_sce)
data(demo_airr)

## -----------------------------------------------------------------------------
vignette("vignette_reproduce_original")

## -----------------------------------------------------------------------------
set.seed(123)

## -----------------------------------------------------------------------------
contig.list <- loadContigs(input = demo_airr, format = "AIRR")

# Format to `scRepertoire`'s requirements and some light filtering
combined.TCR <- combineTCR(contig.list,
    removeNA = TRUE,
    removeMulti = FALSE,
    filterMulti = TRUE
)

## -----------------------------------------------------------------------------
sce <- combineExpression(combined.TCR, demo_sce)

## -----------------------------------------------------------------------------
sce <- setupVdjPseudobulk(sce,
    mode_option = "abT",
    already.productive = TRUE,
    subsetby = "anno_lvl_2_final_clean",
    groups = c("CD8+T", "CD4+T", "ABT(ENTRY)", "DP(P)_T", "DP(Q)_T")
)

## -----------------------------------------------------------------------------
head(colData(sce))

## -----------------------------------------------------------------------------
plotUMAP(sce, color_by = "anno_lvl_2_final_clean")

## ----warning = FALSE----------------------------------------------------------
library(miloR)
milo_object <- Milo(sce)
milo_object <- buildGraph(milo_object, k = 30, d = 20, reduced.dim = "X_scvi")
milo_object <- makeNhoods(milo_object,
    reduced_dims = "X_scvi", d = 20,
    prop = 0.3
)

## ----warning = FALSE----------------------------------------------------------
milo_object <- miloUmap(milo_object, n_neighbors = 30)

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

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

## -----------------------------------------------------------------------------
pb.milo

## -----------------------------------------------------------------------------
pb.milo <- runPCA(pb.milo, assay.type = "Feature_space", ncomponents = 20)
plotPCA(pb.milo, color_by = "anno_lvl_2_final_clean")

## -----------------------------------------------------------------------------
# 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
dm <- DiffusionMap(t(pca), n_pcs = 10, 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",
    knn = 30
)

## -----------------------------------------------------------------------------
head(colData(pb.milo))

## ----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, branch.tips)

## ----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()