## ----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 = 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") gmt_in <- pathwayPCA::read_gmt(gmt_fname) 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) library(scater) # Let's see the original clusters plotReducedDim(mini_pbmc3k, dimred="UMAP", colour_by="clusters") library(gridExtra) p1 <- plotReducedDim(mini_pbmc3k, dimred="UMAP", colour_by="mod_HLA-DRB1") p2 <- plotReducedDim(mini_pbmc3k, dimred="UMAP", colour_by="HLA-DRB1") p3 <- plotReducedDim(mini_pbmc3k, dimred="UMAP", colour_by="mod_CD79A") p4 <- plotReducedDim(mini_pbmc3k, dimred="UMAP", colour_by="CD79A") 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)