## ----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) ## ----Checking module headers, warning=FALSE----------------------------------- mod_names(fc) ## ----Plotting nets------------------------------------------------------------ mod_names(fc) network_plots <- show_net(fc) network_plots[["S100A8"]] network_plots[["MS4A1"]] ## ----------------------------------------------------------------------------- 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 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) ## ---- 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)) } } ## ---- 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) ## ----------------------------------------------------------------------------- TUBB1 <- FeaturePlot(pbmc_small, "TUBB1") DRB1 <- FeaturePlot(pbmc_small, "HLA-DRB1") grid.arrange(TUBB1, DRB1, ncol = 2)