## ----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=6, fig.height=4------------------------------------------------
ggplot(data = d$posterior_summary$epsilon)+
  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(epsilon))+
  scale_x_continuous(expand = c(0,0))

## ----fig.width=4, fig.height=3------------------------------------------------
ggplot(data = d$posterior_summary$delta)+
    facet_wrap(facets = ~contrast, ncol = 1)+
    geom_histogram(aes(mean), bins = 100)+
    xlab(label = expression(delta))+
    theme_bw(base_size = 10)

## ----fig.width=4, fig.height=3------------------------------------------------
ggplot(data = d$posterior_summary$epsilon)+
    facet_wrap(facets = ~contrast, ncol = 1)+
    geom_histogram(aes(mean), bins = 100)+
    xlab(label = expression(epsilon))+
    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))

## ----fig.width=6, fig.height=7------------------------------------------------
epsilon <- d$posterior_summary$epsilon
epsilon <- epsilon[epsilon$contrast %in% c("a-b", "a-c", "b-c"), ]

ggplot(data = epsilon)+
  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(epsilon))+
  scale_x_continuous(expand = c(0,0))

## -----------------------------------------------------------------------------
ggplot(data = d$posterior_summary$delta)+
    facet_wrap(facets = ~contrast, ncol = 1)+
    geom_histogram(aes(mean), bins = 100)+
    xlab(label = expression(delta))+
    theme_bw(base_size = 10)

## -----------------------------------------------------------------------------
ggplot(data = d$posterior_summary$epsilon)+
    facet_wrap(facets = ~contrast, ncol = 1)+
    geom_histogram(aes(mean), bins = 100)+
    xlab(label = expression(epsilon))+
    theme_bw(base_size = 10)

## ----session_info-------------------------------------------------------------
utils::sessionInfo()