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

## ----eval=FALSE---------------------------------------------------------------
# library(devtools)
# install_github("snaketron/ClustIRR", upgrade = "never")

## ----eval = FALSE-------------------------------------------------------------
# library(ClustIRR)

## ----eval = FALSE-------------------------------------------------------------
# library(ggplot2)
# library(ggforce)
# library(igraph)
# library(patchwork)
# library(ggrepel)
# library(rBLAST)

## ----eval = FALSE-------------------------------------------------------------
# # load the data
# data(D_CMV_FLU)

## ----eval = FALSE-------------------------------------------------------------
# # TCR repertoires
# s <- D_CMV_FLU$s
# 
# # Meta data
# meta <- D_CMV_FLU$meta

## ----eval = FALSE-------------------------------------------------------------
# str(s)

## ----eval = FALSE-------------------------------------------------------------
# str(meta)

## ----eval = FALSE-------------------------------------------------------------
# # number of clonotypes
# table(s$sample)
# 
# # number of unique CDR3b sequences
# length(unique(s$CDR3b))
# 
# # number of unique CDR3b sequences
# length(unique(s$CDR3a))
# 
# # count number of cells per samples
# print(aggregate(clone_size~sample, data = s, FUN = sum))

## ----eval = FALSE, fig.width=4, fig.height=2.5--------------------------------
# ggplot(data = s)+
#   geom_sina(aes(x = sample, y = clone_size), shape = 21)+
#   theme_bw()+
#     scale_y_log10()+
#     annotation_logticks(base = 10, sides = "l")

## ----eval = FALSE-------------------------------------------------------------
# # Single cell gene expression data
# gex <- D_CMV_FLU$gex

## ----eval = FALSE-------------------------------------------------------------
# str(gex)

## ----eval = FALSE, eval = FALSE-----------------------------------------------
# o <- clustirr(s = s, meta = meta,
#               control = list(blast_gmi = 0.9)) # min. BLAST sequence identity

## ----eval = FALSE-------------------------------------------------------------
# str(o)

## ----eval = FALSE-------------------------------------------------------------
# is_igraph(o$graph)

## ----eval = FALSE, fig.width=6, fig.height=6----------------------------------
# plot_graph(g = o, as_visnet = T, node_opacity = 0.75, show_singletons = F)

## ----eval = FALSE-------------------------------------------------------------
# n <- as_data_frame(x = o$graph, what = "vertices")
# nrow(n)

## ----eval = FALSE-------------------------------------------------------------
# # get a data.frame 'e' with of edges and their attributes
# e <- as_data_frame(x = o$graph, what = "edges")
# 
# # max edge weight
# e_max <- e[which.max(e$nweight),]
# # get CDR3 sequences
# n$CDR3b[which(n$name %in% c(e_max$from, e_max$to))]

## ----eval = FALSE-------------------------------------------------------------
# # min edge weight
# e_min <- e[which.min(e$nweight),]
# # get CDR3 sequences
# n$CDR3b[which(n$name %in% c(e_min$from, e_min$to))]

## ----eval = FALSE-------------------------------------------------------------
# gcd <- detect_communities(graph = o$graph,
#                           algorithm = "leiden",
#                           chains = "CDR3b")

## ----eval = FALSE-------------------------------------------------------------
# names(gcd)

## ----eval = FALSE-------------------------------------------------------------
# dim(gcd$community_occupancy_matrix)

## ----eval = FALSE-------------------------------------------------------------
# head(gcd$community_occupancy_matrix)

## ----eval = FALSE, fig.width=4.5, fig.height=3--------------------------------
# honeycomb <- get_honeycombs(com = gcd$community_occupancy_matrix)
# honeycomb[[1]]

## ----eval = FALSE-------------------------------------------------------------
# dco_obj <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix,
#                mcmc_control = list(mcmc_warmup = 500,
#                                    mcmc_iter = 1500,
#                                    mcmc_chains = 4,
#                                    mcmc_cores = 4))

## ----eval = FALSE, fig.width=4, fig.height=3----------------------------------
# ggplot(data = dco_obj$posterior_summary$y_hat)+
#   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()

## ----eval = FALSE, fig.width=5, fig.height=3, warning=FALSE-------------------
# bv_cmv <- get_beta_violin_ag(node_summary = gcd$node_summary,
#                              beta = dco_obj$posterior_summary$beta,
#                              ag_species = TRUE,
#                              ag = "CMV",
#                              db = "vdjdb")
# bv_cmv$plot

## ----eval = FALSE, fig.width=5, fig.height=3, warning=FALSE-------------------
# bv_flu <- get_beta_violin_ag(node_summary = gcd$node_summary,
#                              beta = dco_obj$posterior_summary$beta,
#                              ag_species = TRUE,
#                              ag = "InfluenzaA",
#                              db = "vdjdb")
# bv_flu$plot

## ----eval = FALSE-------------------------------------------------------------
# bv_cmv$d[which(bv_cmv$d$sample == "CMV_A0201" &
#                    bv_cmv$d$mean>6.5 &
#                    bv_cmv$d$hit==T),]
# 
# # community 9

## ----eval = FALSE-------------------------------------------------------------
# # Large orange dots with high beta in repertoire Flu
# bv_flu$d[which(bv_flu$d$sample == "Flu" &
#                    bv_flu$d$mean>3 &
#                    bv_flu$d$hit==T),]
# 
# # communities 34, 91, 94

## ----eval = FALSE-------------------------------------------------------------
# ns <- gcd$node_summary
# nrow(ns[ns$community == 9,])

## ----eval = FALSE, fig.width=6, fig.height=1.3--------------------------------
# get_cdr3_motifs(node_summary = gcd$node_summary, community_id = 9)

## ----eval = FALSE, fig.width=5, fig.height=5, warning=FALSE-------------------
# sg <- subgraph(graph = gcd$graph, vids = V(gcd$graph)$community %in% 9)
# 
# V(sg)$Ag <- ifelse(test = regexec(pattern = "CMV", text = V(sg)$Ag_species)!=-1, yes = "+", no = "-")
# V(sg)$col <- ifelse(test = V(sg)$Ag=="+", yes = "orange", no = "#A1BFD9")
# V(sg)$label <- paste0("Community:", V(sg)$community, " (", V(sg)$clone_size, ")" ,
#                       "\n CDR3a:",V(sg)$CDR3a, "\n CDR3b:",V(sg)$CDR3b)
# 
# V(sg)$label[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)] <-
#     paste0("(", V(sg)$clone_size, ")")[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)]
# V(sg)$shape <- ifelse(test = V(sg)$sample=="Flu", yes = "circle", no = "square")
# 
# par(mai=c(0,0,0,0))
# plot(sg,
#      vertex.label = V(sg)$label,
#      layout = igraph::layout_components,
#      vertex.label.cex = 0.8,
#      vertex.size = log(V(sg)$clone_size)+3,
#      vertex.frame.color = NA,
#      vertex.color = V(sg)$col,
#      vertex.shape = V(sg)$shape,
#      edge.width = 0.2)

## ----eval = FALSE-------------------------------------------------------------
# d <- igraph::as_data_frame(x = sg, what = "vertices")
# d[order(d$clone_size, decreasing = T), ][1:2, c("CDR3a", "CDR3b", "clone_size",
#                                                 "sample", "Ag_species")]

## ----eval = FALSE, fig.width=5.5, fig.height=3--------------------------------
# # show the CDR3 logos
# wrap_plots(lapply(X = c(34, 91, 94),
#                   FUN = get_cdr3_motifs, node_summary = gcd$node_summary), ncol = 1)

## ----eval = FALSE, fig.width=10, fig.height=10, warning=FALSE-----------------
# # Show the networks
# sg <- subgraph(graph = gcd$graph, vids = V(gcd$graph)$community %in% c(34, 91, 94))
# 
# V(sg)$Ag <- ifelse(test = regexec(pattern = "CMV", text = V(sg)$Ag_species)!=-1, yes = "+", no = "-")
# V(sg)$col <- ifelse(test = V(sg)$Ag=="+", yes = "orange", no = "#A1BFD9")
# V(sg)$label <- paste0("Community:", V(sg)$community, " (", V(sg)$clone_size, ")" ,
#                       "\n CDR3a:",V(sg)$CDR3a, "\n CDR3b:",V(sg)$CDR3b)
# 
# V(sg)$label[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)] <-
#     paste0("(", V(sg)$clone_size, ")")[which(ave(V(sg)$clone_size, V(sg)$community, FUN = function(x) x == max(x))==0)]
# V(sg)$shape <- ifelse(test = V(sg)$sample=="Flu", yes = "circle", no = "square")
# 
# par(mai=c(0,0,0,0))
# plot(sg,
#      vertex.label = V(sg)$label,
#      layout = igraph::layout_components,
#      vertex.label.cex = 0.8,
#      vertex.size = log(V(sg)$clone_size)+3,
#      vertex.frame.color = NA,
#      vertex.color = V(sg)$col,
#      vertex.shape = V(sg)$shape,
#      edge.width = 0.2)

## ----eval = FALSE-------------------------------------------------------------
# # this is the data.frame with all clonotypes
# ns <- gcd$node_summary
# 
# # extract the clonotypes in community 9 from repertoire CMV_A0201
# ns_com_9 <- ns[ns$sample == "CMV_A0201" & ns$community == 9,]
# 
# # extract the cell_ids (barcodes)
# cell_ids <- unlist(strsplit(x = ns_com_9$cell_id, split = "\\,"))

## ----eval = FALSE-------------------------------------------------------------
# # extract the gex profiles of cells in community 9
# gex_com_9 <- gex[gex$sample == "CMV_A0201" & gex$cell_id %in% cell_ids,]
# # ... and in cells not in community 9
# gex_com_not_9 <- gex[gex$sample == "CMV_A0201" & !gex$cell_id %in% cell_ids,]

## ----eval = FALSE, fig.width=5.5, fig.height=3.5------------------------------
# # compute the sum of gene expressions within each cell
# t_act_9 <- apply(X = gex_com_9[,1:12], MARGIN = 1, FUN = sum)
# t_act_not_9 <- apply(X = gex_com_not_9[,1:12], MARGIN = 1, FUN = sum)
# 
# # compare the distributions
# # red = 9
# # black = not 9
# # do we see a clear difference? yes: red distribution shifted towards larger
# # values -> up-regulation of genes involved in T-cell activation
# plot(density(t_act_9), col = "red",
#      main = "T cell activation in community 9 (red) vs. rest (black)",
#      xlab = "T cell activation score")
# points(density(t_act_not_9), type = "l", col = "black")

## ----eval = FALSE-------------------------------------------------------------
# # this is the data.frame with all clonotypes
# ns <- gcd$node_summary
# 
# # extract the clonotypes in community 34 from repertoire Flu
# ns_com_34 <- ns[ns$sample == "Flu" & ns$community == 34,]
# ns_com_91 <- ns[ns$sample == "Flu" & ns$community == 91,]
# ns_com_94 <- ns[ns$sample == "Flu" & ns$community == 94,]
# 
# # extract the cell_ids (barcodes)
# cell_ids_34 <- unlist(strsplit(x = ns_com_34$cell_id, split = "\\,"))
# cell_ids_91 <- unlist(strsplit(x = ns_com_91$cell_id, split = "\\,"))
# cell_ids_94 <- unlist(strsplit(x = ns_com_94$cell_id, split = "\\,"))

## ----eval = FALSE-------------------------------------------------------------
# # extract the gex profiles of cells in each community
# gex_com_34 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_34,]
# gex_com_91 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_91,]
# gex_com_94 <- gex[gex$sample == "Flu" & gex$cell_id %in% cell_ids_94,]
# 
# # ... and in cells not in communities 34, 91, and 94
# gex_com_not <- gex[gex$sample == "Flu" & !gex$cell_id %in% c(cell_ids_34,
#                                                               cell_ids_91,
#                                                               cell_ids_94),]

## ----eval = FALSE, fig.width=5.4, fig.height=4--------------------------------
# # compute the sum of gene expressions within each cell
# t_act_34 <- apply(X = gex_com_34[,1:12], MARGIN = 1, FUN = sum)
# t_act_91 <- apply(X = gex_com_91[,1:12], MARGIN = 1, FUN = sum)
# t_act_94 <- apply(X = gex_com_94[,1:12], MARGIN = 1, FUN = sum)
# 
# t_act_not <- apply(X = gex_com_not[,1:12], MARGIN = 1, FUN = sum)
# 
# # compare the distributions
# plot(density(t_act_34), col = "pink",
#      main = "T cell activation in community 34, 91, 94, vs. rest",
#      xlab = "T cell activation score")
# points(density(t_act_91), type = "l", col = "forestgreen")
# points(density(t_act_94), type = "l", col = "green")
# points(density(t_act_not), type = "l", col = "black")