---
title: "Finding **biological condition**-specific changes in T- and B-cell 
receptor repertoires with ClustIRR"
output:
    BiocStyle::html_document
vignette: >
    %\VignetteIndexEntry{Finding biological condition-specific changes in T- and B-cell receptor repertoires with ClustIRR}
    %\VignetteEncoding{UTF-8}
    %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r, include=FALSE}
knitr::opts_chunk$set(comment = "", warning = FALSE, message = FALSE)
```

```{r}
library(knitr)
library(ClustIRR)
library(ggplot2)
library(patchwork)
library(ggforce)
library(rBLAST)
theme_set(new = theme_bw(base_size = 10))
```


```{r, include = FALSE}
# only run code if blast+ is installed
run <- has_blast()
```

```{r, echo = FALSE, eval = !run, results='asis'}
cat("**Note: BLAST was not installed when this vignette was built!**")
```


```{r, eval=run}
# 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))
```


```{r, eval=run}
data("D2", package = "ClustIRR")
```

```{r, eval=run}
cl <- clustirr(s = D2, control = list(blast_gmi = 0.8))
```


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


```{r, eval=run}
dim(gcd$community_occupancy_matrix)
```


```{r, fig.width=8, fig.height=9, eval=run}
honeycomb <- get_honeycombs(com = gcd$community_occupancy_matrix)
wrap_plots(honeycomb, nrow = 5, ncol = 3)+
    plot_annotation(tag_levels = 'A')
```



```{r, eval=run}
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))
```



```{r, fig.width=6, fig.height=3, eval=run}
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))
```


```{r, fig.width=7, fig.height=3, eval=run}
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))
```