## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(comment = "", warning = FALSE, message = FALSE)
## -----------------------------------------------------------------------------
library(knitr)
library(ClustIRR)
library(igraph)
library(ggplot2)
library(ggrepel)
library(patchwork)
## ----eval=FALSE---------------------------------------------------------------
# if(!require("BiocManager", quietly = TRUE)) {
# install.packages("BiocManager")
# }
# BiocManager::install("ClustIRR")
## ----graphic, echo = FALSE, fig.align="left", out.width='100%'----------------
knitr::include_graphics("../inst/extdata/logo.png")
## -----------------------------------------------------------------------------
data("CDR3ab", package = "ClustIRR")
## -----------------------------------------------------------------------------
set.seed(127)
n <- 300
# get 300 CDR3a and CDR3b pairs from the data -> IRR a
a <- data.frame(CDR3a = CDR3ab$CDR3a[1:n],
CDR3b = CDR3ab$CDR3b[1:n],
clone_size=rpois(n = n, lambda = 3)+1,
sample = "a")
# get 200 CDR3a and CDR3b pairs from the data -> IRR b
b <- data.frame(CDR3a = CDR3ab$CDR3a[101:(n+100)],
CDR3b = CDR3ab$CDR3b[101:(n+100)],
clone_size=rpois(n = n, lambda = 3)+1,
sample = "b")
## -----------------------------------------------------------------------------
str(a)
## -----------------------------------------------------------------------------
str(b)
## -----------------------------------------------------------------------------
# perform clust_irr analysis of repertoire a
c_a <- cluster_irr(s = a, control = list(gmi = 0.7))
# ... and b
c_b <- cluster_irr(s = b, control = list(gmi = 0.7))
## -----------------------------------------------------------------------------
kable(head(c_a@clust$CDR3a), digits = 2)
## -----------------------------------------------------------------------------
kable(head(c_a@clust$CDR3a), digits = 2)
## -----------------------------------------------------------------------------
# control = list(gmi = 0.7, trim_flank_aa = 3, db_dist = 0, db_custom = NULL)
kable(head(get_clustirr_inputs(c_a)$s), digits = 2)
## -----------------------------------------------------------------------------
g <- get_graph(clust_irr = c_a)
## -----------------------------------------------------------------------------
plot_graph(g, as_visnet = TRUE)
## ----fig.width=6, fig.height=4.5----------------------------------------------
# data.frame of edges and their attributes
e <- igraph::as_data_frame(x = g$graph, what = "edges")
## -----------------------------------------------------------------------------
kable(head(e), digits = 2)
## ----fig.width=5, fig.height=3.5----------------------------------------------
ggplot(data = e)+
geom_density(aes(ncweight, col = chain))+
geom_density(aes(nweight, col = chain), linetype = "dashed")+
theme_bw()+
xlab(label = "edge weight (solid = ncweight, dashed = nweight)")+
theme(legend.position = "top")
## ----fig.width=6, fig.height=6------------------------------------------------
g <- get_joint_graph(clust_irrs = c(c_a, c_b))
## ----fig.width=6, fig.height=6------------------------------------------------
plot_graph(g = g, as_visnet = TRUE, node_opacity = 0.8)
## -----------------------------------------------------------------------------
gcd <- detect_communities(graph = g$graph,
algorithm = "leiden",
resolution = 1,
weight = "ncweight",
metric = "average",
chains = c("CDR3a", "CDR3b"))
## -----------------------------------------------------------------------------
names(gcd)
## -----------------------------------------------------------------------------
dim(gcd$community_occupancy_matrix)
## -----------------------------------------------------------------------------
head(gcd$community_occupancy_matrix)
## ----fig.width=5, fig.height=5------------------------------------------------
plot(x = gcd$community_occupancy_matrix[, "a"],
y = gcd$community_occupancy_matrix[, "b"],
xlab = "community occupancy [cells in IRR a]",
ylab = "community occupancy [cells in IRR b]")
## -----------------------------------------------------------------------------
kable(head(gcd$community_summary), digits = 2)
## ----fig.width=6, fig.height=4------------------------------------------------
ggplot(data = gcd$community_summary)+
geom_point(aes(x = w_CDR3a, y = w_CDR3b, size = cells_n), shape=21)+
xlab(label = "CDR3a ncweight")+
ylab(label = "CDR3b ncweight")+
scale_size_continuous(range = c(0.5, 5))+
theme_bw(base_size = 10)
## -----------------------------------------------------------------------------
kable(head(gcd$node_summary), digits = 2)
## -----------------------------------------------------------------------------
d <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix,
mcmc_control = list(mcmc_warmup = 500,
mcmc_iter = 1500,
mcmc_chains = 3,
mcmc_cores = 1,
mcmc_algorithm = "NUTS",
adapt_delta = 0.9,
max_treedepth = 10))
## ----fig.width=6, fig.height=2.5----------------------------------------------
ggplot(data = d$posterior_summary$y_hat)+
facet_wrap(facets = ~sample, nrow = 1, scales = "free")+
geom_abline(intercept = 0, slope = 1, linetype = "dashed", col = "gray")+
geom_errorbar(aes(x = y_obs, y = mean, ymin = L95, ymax = H95),
col = "darkgray", width=0)+
geom_point(aes(x = y_obs, y = mean), size = 0.8)+
xlab(label = "observed y")+
ylab(label = "predicted y (and 95% HDI)")+
theme_bw(base_size = 10)
## ----fig.width=6, fig.height=4------------------------------------------------
ggplot(data = d$posterior_summary$delta)+
facet_wrap(facets = ~contrast, ncol = 1)+
geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95),
col = "darkgray", width =0)+
geom_point(aes(x = community, y = mean), size = 0.5)+
theme_bw(base_size = 10)+
theme(legend.position = "none")+
ylab(label = expression(delta))+
scale_x_continuous(expand = c(0,0))
## ----fig.width=4, fig.height=3------------------------------------------------
ggplot(data = d$posterior_summary$delta)+
geom_histogram(aes(mean), bins = 100)+
xlab(label = expression(delta))+
theme_bw(base_size = 10)
## ----echo=FALSE, include=FALSE------------------------------------------------
rm(a, b, c_a, c_b, d, e, g, gcd, n)
## -----------------------------------------------------------------------------
# repertoire size
n <- 200
# a
a <- data.frame(CDR3a = CDR3ab$CDR3a[1:n],
CDR3b = CDR3ab$CDR3b[1:n],
sample = "a")
# b
b <- data.frame(CDR3a = CDR3ab$CDR3a[1:n],
CDR3b = CDR3ab$CDR3b[1:n],
sample = "b")
# c
c <- data.frame(CDR3a = CDR3ab$CDR3a[1:n],
CDR3b = CDR3ab$CDR3b[1:n],
sample = "c")
## -----------------------------------------------------------------------------
get_clonal_expansion <- function(n, p_expanded) {
s <- sample(x = c(0, 1), size = n, prob = c(1-p_expanded,
p_expanded), replace = T)
y <- vapply(X = s, FUN.VALUE = numeric(1), FUN = function(x) {
if(x == 0) {
return(rpois(n = 1, lambda = 0.5))
}
return(rpois(n = 1, lambda = 50))
})
return(y)
}
## -----------------------------------------------------------------------------
# simulate expansion of specific communities
set.seed(1243)
clone_size <- rpois(n = n, lambda = 3)+1
expansion_factor <- get_clonal_expansion(n = n, p_expanded = 0.02)
a$clone_size <- clone_size
b$clone_size <- clone_size+expansion_factor*1
c$clone_size <- clone_size+expansion_factor*2
## -----------------------------------------------------------------------------
# run cluster_irr on each IRR and join the results
clust_irrs <- c(cluster_irr(s = a), cluster_irr(s = b), cluster_irr(s = c))
## -----------------------------------------------------------------------------
g <- get_joint_graph(clust_irrs = clust_irrs, cores = 1)
## ----fig.width=6, fig.height=6------------------------------------------------
plot_graph(g = g, as_visnet = TRUE, node_opacity = 0.8)
## -----------------------------------------------------------------------------
gcd <- detect_communities(graph = g$graph,
weight = "ncweight",
algorithm = "leiden",
resolution = 1,
chains = c("CDR3a", "CDR3b"))
## ----fig.width=7, fig.height=2.5----------------------------------------------
(ggplot(data = gcd$community_summary)+
geom_point(aes(x = cells_a, y = cells_b), shape = 21)+
theme_bw(base_size = 10))|
(ggplot(data = gcd$community_summary)+
geom_point(aes(x = cells_a, y = cells_c), shape = 21)+
theme_bw(base_size = 10))|
(ggplot(data = gcd$community_summary)+
geom_point(aes(x = cells_b, y = cells_c), shape = 21)+
theme_bw(base_size = 10))+
plot_annotation(tag_level = 'A')
## -----------------------------------------------------------------------------
community_occupancy_matrix <- gcd$community_occupancy_matrix
head(community_occupancy_matrix)
## -----------------------------------------------------------------------------
d <- dco(community_occupancy_matrix = community_occupancy_matrix,
mcmc_control = list(mcmc_warmup = 300,
mcmc_iter = 900,
mcmc_chains = 3,
mcmc_cores = 1,
mcmc_algorithm = "NUTS",
adapt_delta = 0.95,
max_treedepth = 11))
## ----fig.width=6, fig.height=2.5----------------------------------------------
ggplot(data = d$posterior_summary$y_hat)+
facet_wrap(facets = ~sample, nrow = 1, scales = "free")+
geom_abline(intercept = 0, slope = 1, linetype = "dashed", col = "gray")+
geom_errorbar(aes(x = y_obs, y = mean, ymin = L95, ymax = H95),
col = "darkgray", width=0)+
geom_point(aes(x = y_obs, y = mean), size = 0.8)+
xlab(label = "observed y")+
ylab(label = "predicted y (and 95% HDI)")+
theme_bw(base_size = 10)
## ----fig.width=6, fig.height=7------------------------------------------------
beta <- d$posterior_summary$beta
ggplot(data = beta)+
facet_wrap(facets = ~sample, ncol = 1)+
geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95,
col = sample), width = 0)+
geom_point(aes(x = community, y = mean, col = sample), size = 0.5)+
theme_bw(base_size = 10)+
theme(legend.position = "top")+
ylab(label = expression(beta))+
scale_x_continuous(expand = c(0,0))
## ----fig.width=6, fig.height=7------------------------------------------------
delta <- d$posterior_summary$delta
delta <- delta[delta$contrast %in% c("a-b", "a-c", "b-c"), ]
ggplot(data = delta)+
facet_wrap(facets = ~contrast, ncol = 1)+
geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95), width=0)+
geom_point(aes(x = community, y = mean), size = 0.5)+
theme_bw(base_size = 10)+
theme(legend.position = "top")+
ylab(label = expression(delta))+
scale_x_continuous(expand = c(0,0))
## ----session_info-------------------------------------------------------------
utils::sessionInfo()