---
title: 'fcoex: co-expression for single-cell data integrated with Seurat'
author:
- name: Tiago Lubiana
  affiliation: Computational Systems Biology Laboratory, University of São Paulo, Brazil
output:
  html_document:
    toc: yes
  pdf_document:
    toc: yes
  prettydoc::html_pretty:
    highlight: github
    theme: cayman
package: fcoex
vignette: > 
  %\VignetteIndexEntry{fcoex: co-expression for single-cell data integrated with Seurat}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---
## Integration to Seurat 

This vignette is targeted at Seurat users that want to use fcoex. 

A more in depth explanation of each fcoex step is available in the general vignette, which is based on bioconductor packages.

To integrate a pre-processed Seurat object with fcoex, it is only necessary to obtain the normalized expression matrix and the cluster identities, and then run the pipeline:

```{r Running Seurat pipeline, warning=FALSE}
library(Seurat)
library(fcoex)
library(ggplot2)

data(pbmc_small)

exprs <- data.frame(GetAssayData(pbmc_small))
target <- Idents(pbmc_small)

fc <- new_fcoex(data.frame(exprs),target)
fc <- discretize(fc)
fc <- find_cbf_modules(fc,n_genes = 70, verbose = FALSE, is_parallel = FALSE)
fc <- get_nets(fc)
```
Let's take a look at the modules headers:

```{r Checking module headers, warning=FALSE}
mod_names(fc)
```

S100A8 is a marker of some granulocytes and MS4A1 is a markers of B cells. Let's take a look at those two modules. 

```{r Plotting nets}
mod_names(fc)

network_plots <- show_net(fc)

network_plots[["S100A8"]]
network_plots[["MS4A1"]]
```


Once again we can use the Reactome pathways from the package CEMiTool to exemplify how to run an enrichment for human genes.
It is likely that some modules will not have any enrichment, leading to messages of the type "no gene can be mapped.". 
That is not a problem.
If you are not working with human genes, you can just skip this part.

```{r}
gmt_fname <- system.file("extdata", "pathways.gmt", package = "CEMiTool")
gmt_in <- pathwayPCA::read_gmt(gmt_fname)
fc <- mod_ora(fc, gmt_in)

# In Seurat's sample data, pbmc small, no enrichments are found. 
# That is way plot_ora is commented out.

# fc <- plot_ora(fc)
```

```{r Saving Seurat plots, eval = FALSE}
save_plots(name = "fcoex_vignette_Seurat", fc, force = TRUE, directory = "./Plots")
```

### Plotting clusters with  Seurat

One application of fcoex is to find overlapping populations of cells in the dataset, by employing a module-based view of the gene expression landscape. 

Why reclustering? Well, the algorithm behind fcoex considers at the same time inverse and direct correlations. 
Thus, it is nonsense to obtain a plot of the average expression of the genes in a module, for example.  
The reclustering allows us to gather all the signal in the dataset, positive and negative, to see how cells in the dataset behave in relation to that subpart of the transcriptomic space.

Now let's recluster the cells for some of the fcoex module and visualize it in Seurat .

The cells will be divided on two groups: header positive (HP) and header negative (HN). 
For details on why, see the next session on anti-correlated genes.

```{r Plotting reclusters}
library(gridExtra)

fc <- recluster(fc) 

# Running UMAP to obtain layout for cells
pbmc_small <- RunUMAP(pbmc_small, dims = 1:10)

plot1 <- DimPlot(pbmc_small)

pbmc_reclustered <- pbmc_small

# S100A8 is a marker of some monocytes

Idents(pbmc_reclustered) <- idents(fc)[["S100A8"]]
plot2  <- DimPlot(pbmc_reclustered, reduction = 'umap', cols = c("darkgreen", "dodgerblue3")) +
  ggtitle("S100A8")
    

# HLA-DPB1 is a marker of Antigen Presenting Cells
Idents(pbmc_reclustered) <- idents(fc)[["HLA-DPB1"]]
plot3  <- DimPlot(pbmc_reclustered, reduction = 'umap', cols = c("darkgreen", "dodgerblue3")) +
  ggtitle("HLA-DPB1")

# MS4A1 is a marker of B cells
Idents(pbmc_reclustered) <- idents(fc)[["MS4A1"]]
plot4  <- DimPlot(pbmc_reclustered, reduction = 'umap', cols = c("darkgreen", "dodgerblue3")) +
  ggtitle("MS4A1")

grid.arrange(plot1, plot2, plot3, plot4, nrow=2)

```
If you notice the "HP" clusters, you will see that the "HLA-DPB1" recluster covers roughly the same cells as the combination of the "S100A8" recluster and the "MS4A1" reclusters. 
This is expected, as both B cells and monocytes are antigen presenting cells. 

Notice that these overlapped insights are not possible looking by one single flat output.

### Detecting anticorrelated genes in the modules

The  fcoex modules capture also negative correlations.
hese can be especially interesting, as they point to complementary cell classes.

Each module has one "header" gene. The other genes in the module are either in the same cells of the header (thus, header-positive cells, or HP), or present in the other cells of the dataset (header negative). 

For all the modules, we will check to see if such anti-correlations are present .

```{r, message=FALSE, warning=FALSE}

for (i in names(module_genes(fc))) {
  Idents(pbmc_reclustered) <-   fc@mod_idents[[i]]
  
  print(paste("Checking for anticorrelation in module", i))
 
  # Identify markers only for module genes:
  module_genes_in_clusters <-
    FindAllMarkers(
      pbmc_reclustered,
      logfc.threshold = 1,
      only.pos = TRUE,
      features = fc@module_list[[i]],
      verbose = FALSE
    )
  
  # If there are markers of the HN cluster, it means that we have anticorrelation
  if ("HN" %in% module_genes_in_clusters$cluster) {
    module_genes_in_clusters$module = i
    message(paste0("anticorrelated genes found for module ", i))
  }
}
```

There seems that only in the module HLA-DRB1 we have negative correlations. Let's visualize then


```{r, message=FALSE, warning=FALSE}
Idents(pbmc_reclustered) <-   fc@mod_idents[["HLA-DRB1"]]

markers_for_cluster_HLA_DRB1 <-
    FindAllMarkers(
      pbmc_reclustered,
      logfc.threshold = 1,
      only.pos = TRUE,
      features = fc@module_list[["HLA-DRB1"]],
      verbose = FALSE
    )

print(markers_for_cluster_HLA_DRB1)
```

We can see that the HLA genes are in the "HP" direction and TUBB1, GP9 and NGFRAP1 are in the opposing direction.

Let's look at the expression patterns in the UMAP plot.

```{r}
TUBB1 <- FeaturePlot(pbmc_small, "TUBB1")
DRB1 <-  FeaturePlot(pbmc_small, "HLA-DRB1")

grid.arrange(TUBB1, DRB1, ncol = 2)

```

The completeness is not perfect (it is a small example dataset), but we can already see the disjoint pattern of expression. 

The anticorrelations provide an additional layer of exploration of the dataset, and might indicate complementary populations of cells, like senders and receivers, or disjoint expression programs.
