## ----Loading datasets, message=FALSE------------------------------------------
library(fcoex, quietly = TRUE)
library(SingleCellExperiment, quietly = TRUE)
data("mini_pbmc3k")


## ----Plotting single-cell object----------------------------------------------
mini_pbmc3k

## ----Creating fcoex object, message=FALSE-------------------------------------
target <- colData(mini_pbmc3k)
target <- target$clusters
exprs <- as.data.frame(assay(mini_pbmc3k, 'logcounts'))

fc <- new_fcoex(data.frame(exprs),target)


## ----Discretizing dataset, message=FALSE--------------------------------------

fc <- discretize(fc, number_of_bins = 8)

## ----Finding cbf modules, message=FALSE---------------------------------------
fc <- find_cbf_modules(fc,n_genes_selected_in_first_step = 200, verbose = FALSE, is_parallel = FALSE)

## ----Plotting module networks, message=FALSE----------------------------------
fc <- get_nets(fc)

# Taking a look at the first two networks: 
show_net(fc)[["CD79A"]]
show_net(fc)[["HLA-DRB1"]]

## ----Saving plots, eval= FALSE, message=FALSE, results='hide'-----------------
#  save_plots(name = "fcoex_vignette", fc,force = TRUE, directory = "./Plots")

## ----Running ORA analysis, warning=FALSE--------------------------------------
gmt_fname <- system.file("extdata", "pathways.gmt", package = "CEMiTool")

if (gmt_fname != "")
  {
    gmt_in <- pathwayPCA::read_gmt(gmt_fname,  description = TRUE)

} else {
    print("You likely need to install CEMiTool")

}
fc <- mod_ora(fc, gmt_in)
fc <- plot_ora(fc)

## ----Saving plots again,  eval= FALSE, message=FALSE, results='hide'----------
#  save_plots(name = "fcoex_vignette", fc, force = TRUE, directory = "./Plots")

## ----Reclustering , message=FALSE---------------------------------------------

fc <- recluster(fc)


## ----Visualizing--------------------------------------------------------------

colData(mini_pbmc3k) <- cbind(colData(mini_pbmc3k), `mod_HLA_DRB1` = idents(fc)$`HLA-DRB1`)
colData(mini_pbmc3k) <- cbind(colData(mini_pbmc3k), mod_CD79A = idents(fc)$CD79A)


# Let's see the original clusters
library(schex)

mini_pbmc3k <- make_hexbin(mini_pbmc3k, nbins = 40, 
    dimension_reduction = "UMAP", use_dims=c(1,2))

plot_hexbin_meta(mini_pbmc3k, col="clusters", action="majority")

library(gridExtra)
p1 = plot_hexbin_feature_plus(mini_pbmc3k,
    col="clusters", type="logcounts",
    feature="CD79A", action="mean") +
  ggtitle("original clusters (CD79A expression)") +
  theme_void()

p2 =plot_hexbin_feature_plus(mini_pbmc3k,
    col="clusters", type="logcounts",
    feature="HLA-DRB1", action="mean") +
  ggtitle("original clusters (HLA-DRB1 expression)") +
  theme_void()

p3 =  plot_hexbin_feature_plus(mini_pbmc3k,
    col="mod_CD79A", type="logcounts",
    feature="CD79A", action="mean") +
  ggtitle("fcoex CD79A clusters (CD79A expression)") +
  theme_void()

p4 =  plot_hexbin_feature_plus(mini_pbmc3k,
    col="mod_HLA_DRB1", type="logcounts",
    feature="HLA-DRB1", action="mean")+
  ggtitle("fcoex HLA cluster (HLA-DRB1 expression)") +
  theme_void()

grid.arrange(p1, p2, p3, p4, nrow=2)

## ----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)

## -----------------------------------------------------------------------------
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)

## ----Saving Seurat plots, eval = FALSE----------------------------------------
#  save_plots(name = "fcoex_vignette_Seurat", fc, force = TRUE, directory = "./Plots")

## ----Plotting and saving reclusters,  eval = FALSE----------------------------
#  
#  fc <- recluster(fc)
#  
#  file_name <- "pbmc3k_recluster_plots.pdf"
#  directory <- "./Plots/"
#  
#  pbmc_small <- RunUMAP(pbmc_small, dims = 1:10)
#  
#  pdf(paste0(directory,file_name), width = 3, height = 3)
#  print(DimPlot(pbmc_small))
#  for (i in names(module_genes(fc))){
#    Idents(pbmc_small ) <-   idents(fc)[[i]]
#    mod_name <- paste0("M", which(names(idents(fc)) == i), " (", i,")")
#  
#    plot2 <- DimPlot(pbmc_small, reduction = 'umap', cols = c("darkgreen", "dodgerblue3")) +
#      ggtitle(mod_name)
#      print(plot2)
#  }
#  dev.off()

## -----------------------------------------------------------------------------
  fc <- recluster(fc) 
  pbmc_small <- RunUMAP(pbmc_small, dims = 1:10)
  
  
  
  Idents(pbmc_small ) <-   target
  p1 <- DimPlot(pbmc_small)
  Idents(pbmc_small ) <-   idents(fc)[["HLA-DRB1"]]
  
  mod_name <- paste0("M", which(names(idents(fc)) == "HLA-DRB1"), " (", "HLA-DRB1",")")

  p2 <- DimPlot(pbmc_small, cols = c("darkgreen", "dodgerblue3")) +
    ggtitle(mod_name) 
  
  # CD79A is a marker of B cells
  CD79A <- FeaturePlot(pbmc_small, "CD79A")

  # AIF1 is a marker of monocytes
  AIF1 <- FeaturePlot(pbmc_small, "AIF1")
  
  
  library(gridExtra)
  grid.arrange(p1, p2, p3,p4, ncol = 2)
  

## ---- message=FALSE-----------------------------------------------------------

for (i in names(module_genes(fc))){
Idents(pbmc_small ) <-   fc@mod_idents[[i]]

# This bit prints which gene in the module belongs to each cluster. 
# HP is the header-positive cluster (containing SOX19A), HN is the header negative cluster (not containing SOX19A)
# The "features = fc@module_list[[i]]" parameter tells Seurat to compare only the genes in the module "i"
# By removing this parameter, you can potentially expand the list that was retrieved originally by fcoex

# Run only for module genes:
module_genes_in_clusters <- FindAllMarkers(pbmc_small, logfc.threshold = 1, only.pos = TRUE, features = fc@module_list[[i]] )

if("HN" %in% module_genes_in_clusters$cluster){
module_genes_in_clusters$module = i
message(paste0("anticorrelated genes found for module ", i))
print(module_genes_in_clusters) 
}
}

## -----------------------------------------------------------------------------
 TUBB1 <- FeaturePlot(pbmc_small, "TUBB1")
 DRB1 <-  FeaturePlot(pbmc_small, "HLA-DRB1")
  
  
  library(gridExtra)
  grid.arrange(p1, p2, TUBB1, DRB1, ncol = 2)