library(knitr)
library(ClustIRR)
library(ggplot2)
library(patchwork)
library(ggforce)
library(rBLAST)
theme_set(new = theme_bw(base_size = 10))

Note: BLAST was not installed when this vignette was built!

# System requirements
#  The BLAST+ software needs to be installed on your system. Installation
# instructions are available in this package's
# [INSTALL](https://github.com/mhahsler/rBLAST/blob/devel/INSTALL) file and
# at \url{https://www.ncbi.nlm.nih.gov/books/NBK569861/}.

# R needs to be able to find the executable. After installing the software,
# try in R
Sys.which("blastp")

# If the command returns "" instead of the path to the executable,
# then you need to set the environment variable called PATH. In R
Sys.setenv(PATH = paste(Sys.getenv("PATH"),
   "path_to_your_BLAST_installation", sep=.Platform$path.sep))

data("D2", package = "ClustIRR")

cl <- clustirr(s = D2, control = list(blast_gmi = 0.8))

gcd <- detect_communities(graph = cl$graph, 
                          algorithm = "leiden",
                          metric = "average",
                          resolution = 1,
                          iterations = 100,
                          chains = c("CDR3a", "CDR3b"))

dim(gcd$community_occupancy_matrix)

honeycomb <- get_honeycombs(com = gcd$community_occupancy_matrix)
wrap_plots(honeycomb, nrow = 5, ncol = 3)+
    plot_annotation(tag_levels = 'A')

d <- dco(community_occupancy_matrix = gcd$community_occupancy_matrix,
         groups = c(1, 1, 1, 2, 2, 2),
         mcmc_control = list(mcmc_warmup = 300,
                             mcmc_iter = 600,
                             mcmc_chains = 2,
                             mcmc_cores = 1,
                             mcmc_algorithm = "NUTS",
                             adapt_delta = 0.9,
                             max_treedepth = 10))

ggplot(data = d$posterior_summary$beta)+
    geom_sina(aes(x = sample, y = mean))|
    ggplot(data = d$posterior_summary$beta_mu)+
    geom_sina(aes(x = as.character(g), y = mean))

ggplot(data = d$posterior_summary$epsilon)+
    facet_wrap(facets = ~contrast, ncol = 2)+
    geom_errorbar(aes(x = community, y = mean, ymin = L95, ymax = H95), 
                  col = "lightgray", width = 0)+
    geom_point(aes(x = community, y = mean), shape = 21, fill = "black", 
               stroke = 0.4, col = "white", size = 1.25)+
    theme(legend.position = "top")+
    ylab(label = expression(delta))+
    scale_x_continuous(expand = c(0,0))