fcoex: co-expression for single-cell data

Computational Systems Biology Laboratory, University of São Paulo, Brazil

Introduction and basic pipeline

The goal of fcoex is to provide a simple and intuitive way to generate co-expression nets and modules for single cell data. It is based in 3 steps:

• Pre-processing and label assignement (prior to fcoex)
• Discretization of gene expression
• Correlation and module detection via the FCBF algorithm (Fast Correlation-Based Filter)

First of all, we will a already preprocessed single cell dataset from 10XGenomics ( preprocessed according to the OSCA pipeline, 14/08/2019). It contains peripheral blood mononuclear cells and the most variable genes.

library(fcoex, quietly = TRUE)

## Warning: replacing previous import 'Matrix::head' by 'utils::head' when loading
## 'fcoex'

library(SingleCellExperiment, quietly = TRUE)
data("mini_pbmc3k")

This is the single cell object we will explore in this vignette:

mini_pbmc3k

## class: SingleCellExperiment 
## dim: 1700 600 
## metadata(0):
## assays(2): counts logcounts
## rownames(1700): LYZ S100A9 ... ALG10 LLGL1
## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
## colnames(600): V1 V2 ... V599 V600
## colData names(14): Sample Barcode ... mod_FCER1G mod_CD3D
## reducedDimNames(2): PCA UMAP
## mainExpName: NULL
## altExpNames(0):

Creating the fcoex object

Now let’s use the normalized data and the cluster labels to build the co-expresison networks. The labels were obtained by louvain clustering on a graph build from nearest neighbours. That means that these labels are a posteriori, and this depends on the choice of the analyst.

The fcoex object is created from 2 different pieces: a previously normalized expression table (genes in rows) and a target factor with classes for the cells.

target <- colData(mini_pbmc3k)
target <- target$clusters
exprs <- as.data.frame(assay(mini_pbmc3k, 'logcounts'))

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

The first step is the conversion of the count matrix into a discretized dataframe. The standar of fcoex is a simple binarization that works as follows:

For each gene, the maximum and minimum values are stored. This range is divided in n bins of equal width (parameter to be set). The first bin is assigned to the class “low” and all the others to the class “high”.

fc <- discretize(fc, number_of_bins = 8)

Getting the modules

Note that many other discretizations are avaible, from the implementations in the FCBF Bioconductor package. This step affects the final results in many ways. However, we found empirically that the default parameter often yields interesting results.

After the discretization, we proceed to constructing a network and extracting modules. The co-expression adjacency matriz generated is modular in its inception. All correlations used are calculated via Symmetrical Uncertainty. Three steps are present:

1 - Selection of n genes to be considered, ranked by correlation to the target variable.

2 - Detection of predominantly correlated genes, a feature selection approach defined in the FCBF algorithm

3 - Building of modules around selected genes. Correlations between two genes are kept if they are more correlated to each other than to the target lables

You can choose either to have a non-parallel processing, with a progress bar, or a faster parallel processing without progress bar. Up to you.

fc <- find_cbf_modules(fc,n_genes_selected_in_first_step = 200, verbose = FALSE, is_parallel = FALSE)

## [1] 0.1
## [1] "Number of features features =  1700"
## [1] "Number of prospective features =  199"
## [1] "Number of final features =  28"
## [1] "FCER1G"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TYROBP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD3D"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DRB1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LGALS2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCN1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IL32"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "AIF1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DRB5"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CFD"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LST1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CST3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IFITM3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A9"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A8"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TYMP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LGALS1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SERPINA1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DPA1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IFI30"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DRA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD79A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CST7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LTB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MS4A6A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD68"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SPI1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MS4A1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DQA1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NKG7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GZMA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LGALS3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD79B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RP11-290F20-3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CTSW"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CFP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD3E"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DPB1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD14"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GZMB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GRN"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A11"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CPVL"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CDA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "BRI3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CTSS"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A4"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TCL1A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCGR3A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CEBPD"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FGFBP2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DMA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "APOBEC3A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CCL5"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PTPRCAP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LY86"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PILRA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DQB1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PRF1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CSF3R"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NCF2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FAM26F"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MS4A7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LINC00926"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PYCARD"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LYZ"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "COTL1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GPX1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GZMH"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LDHB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SAT1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PSAP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SRGN"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "BLVRB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NPC2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GABARAP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GNLY"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TKT"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GSTP1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TIMP1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CEBPB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PLBD1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IL7R"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MAFB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DQA2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TNFSF13B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ANXA2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "VPREB3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ALDH2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GZMM"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "C1orf162"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "AP2S1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RGS2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DMB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "STXBP2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD247"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "AP1S2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RAC1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HOPX"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ASGR1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RAB32"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CCL4"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GPBAR1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A6"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FPR1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCGRT"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "XCL2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RHOC"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "BANK1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CCL3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CTSZ"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IGFBP7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "STX11"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ANXA1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A12"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PRKCQ-AS1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LYN"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PRELID1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MEF2C"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GZMK"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "VSTM1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NEAT1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IFI6"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD8A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCER2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RBP7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MPEG1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NUP214"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD74"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LILRA3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCRLA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "KLRD1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IL1B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ISG15"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GIMAP5"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD63"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HVCN1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SPON2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "VIM"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MATK"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GIMAP7"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HCST"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CLIC3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD69"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "NOSIP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FOLR3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "AKR1C3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD19"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MT2A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "S100A10"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD3G"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TSPAN13"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ID2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IL8"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "APMAP"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RHOB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SCPEP1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RNASE6"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RGCC"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "BLK"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GPR56"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD27"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PRSS23"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MT-CO2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD72"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HLA-DOB"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD8B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "YBX3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LINC00936"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TTC38"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "KIAA0125"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SOD2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "PDLIM1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CD99"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "H2AFY"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "IFITM1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "G0S2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CXCR4"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "LAMTOR2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "STK17A"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SH3BGRL3"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SPOCK2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCGR2B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TNFSF10"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TNFRSF4"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CHST12"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "GAPT"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ZFP36L2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "TIGIT"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "XCL1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "HAVCR2"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "RETN"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "CYBA"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "KLRG1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "ARPC1B"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "MARCH1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "SUPT4H1"

## Warning in xtfrm.data.frame(x): cannot xtfrm data frames

## [1] "FCER1G"
## [1] "CD3D"
## [1] "HLA-DRB1"
## [1] "AIF1"
## [1] "CST3"
## [1] "CD79A"
## [1] "CST7"
## [1] "CTSS"
## [1] "FCGR3A"
## [1] "CD7"
## [1] "COTL1"
## [1] "SAT1"
## [1] "GZMM"
## [1] "CD247"
## [1] "S100A6"
## [1] "ANXA1"
## [1] "NEAT1"
## [1] "VIM"
## [1] "CD69"
## [1] "NOSIP"
## [1] "RGCC"
## [1] "MT-CO2"
## [1] "CD99"
## [1] "IFITM1"
## [1] "CXCR4"
## [1] "STK17A"
## [1] "TNFRSF4"
## [1] "ZFP36L2"

The get_nets function takes the modules and plot networks. The visualizations it generates were heavily inspired by the CEMiTool package, as much of the code in fcoex was.

We will take a look at the first two networks

fc <- get_nets(fc)

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

show_net(fc)[["HLA-DRB1"]]

To save the plots, you can run the save plots function, which will create a “./Plots” directory and store plots there.

save_plots(name = "fcoex_vignette", fc,force = TRUE, directory = "./Plots")

Running an enrichment analysis

You can also run an over-representation analysis to see if the modules correspond to any known biological pathway. In this example we will use the Reactome groups available in the CEMiTool package.

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

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)

## --> No gene can be mapped....

## --> Expected input gene ID: ABCA6,GNB2,MAPK14,SRC,ABCD1,PSENEN

## --> return NULL...

## --> No gene can be mapped....

## --> Expected input gene ID: PSENEN,RPL37,ABCD3,NCAN,RPL9,ABCC6

## --> return NULL...

fc <- plot_ora(fc)

Now we can save the plots again. Note that we have to set the force parameter equal to TRUE now, as the “./Plots” directory was already created in the previous step.

save_plots(name = "fcoex_vignette", fc, force = TRUE, directory = "./Plots")

Reclustering the cells to find module-based populations.

There is now a folder with the correlations, to explore the data.

We will use the module assignments to subdivide the cells in populations of interest. This is a way to explore the data and look for possible novel groupings ignored in the previous clustering step.

fc <- recluster(fc)

## [1] "FCER1G"
## [1] "CD3D"
## [1] "HLA-DRB1"
## [1] "AIF1"
## [1] "CST3"
## [1] "CD79A"
## [1] "CST7"
## [1] "CTSS"
## [1] "FCGR3A"
## [1] "CD7"
## [1] "SAT1"
## [1] "S100A6"

We generated new labels based on each fcoex module. Now we will visualize them using UMAP. Let’s see the population represented in the modules CD79A and HLA-DRB1. Notably, the clustering patterns are largely influenced by the expression patterns of header genes. It is interesting to see that two groups are present, header-positive (HP) and header negative (HN) clusters.

The stratification and exploration of different clustering points of view is one of the core features of fcoex.

Plotting clusters with scater

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)

## Loading required package: Seurat

## Attaching SeuratObject

## 
## Attaching package: 'Seurat'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     Assays

## Loading required package: ggplot2

## Loading required package: shiny

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)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:Biobase':
## 
##     combine

## The following object is masked from 'package:BiocGenerics':
## 
##     combine

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)

Integration to Seurat

Another popular framework for single-cell analysis is Seurat. Let’s see how to integrate it with fcoex.

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)

## Created new fcoex object.

fc <- discretize(fc)
fc <- find_cbf_modules(fc,n_genes = 70, verbose = FALSE, is_parallel = FALSE)

## [1] 0.1

## Getting SU scores

## Running FCBF to find module headers

## [1] "Number of features features =  230"
## [1] "Number of prospective features =  69"
## [1] "Number of final features =  9"

## Calculating adjacency matrix

## [1] "S100A8"
## [1] "HLA-DPB1"
## [1] "TYMP"
## [1] "HLA-DRB1"
## [1] "HLA-DPA1"
## [1] "AIF1"
## [1] "LST1"
## [1] "HLA-DRB5"
## [1] "FCGRT"
## [1] "HLA-DQB1"
## [1] "IFITM3"
## [1] "CST3"
## [1] "S100A9"
## [1] "S100A11"
## [1] "FCER1G"
## [1] "TYROBP"
## [1] "FCN1"
## [1] "CD14"
## [1] "SMCO4"
## [1] "SERPINA1"
## [1] "GRN"
## [1] "HLA-DQA1"
## [1] "IFI30"
## [1] "MS4A6A"
## [1] "MS4A1"
## [1] "PSAP"
## [1] "HLA-DMB"
## [1] "FPR1"
## [1] "LGALS3"
## [1] "LGALS1"
## [1] "BID"
## [1] "LYZ"
## [1] "C5AR1"
## [1] "HCK"
## [1] "CFD"
## [1] "CTSB"
## [1] "CD7"
## [1] "IFI6"
## [1] "HLA-DQA2"
## [1] "CTSS"
## [1] "CFP"
## [1] "CD79A"
## [1] "NUP214"
## [1] "LGALS2"
## [1] "HLA-DRA"
## [1] "HLA-DMA"
## [1] "SAT1"
## [1] "NFKBIA"
## [1] "CARD16"
## [1] "TCL1A"
## [1] "IL1B"
## [1] "RNF130"
## [1] "LY86"
## [1] "ASGR1"
## [1] "LINC00926"
## [1] "FCER2"
## [1] "GZMM"
## [1] "LINC00936"
## [1] "TSPO"
## [1] "GSTP1"
## [1] "HVCN1"
## [1] "GP9"
## [1] "TUBB1"
## [1] "NGFRAP1"
## [1] "LRRC25"
## [1] "LCK"
## [1] "TALDO1"
## [1] "GZMA"
## [1] "GNLY"

## Trimming and getting modules from adjacency matrix

## [1] "S100A8"
## [1] "HLA-DPB1"
## [1] "TYMP"
## [1] "HLA-DRB1"
## [1] "MS4A1"
## [1] "LYZ"
## [1] "CD7"
## [1] "IFI6"
## [1] "HLA-DQA2"

fc <- get_nets(fc)

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.

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)

save_plots(name = "fcoex_vignette_Seurat", fc, force = TRUE, directory = "./Plots")

Plotting clusters with Seurat

Now let’s recluster fcoex and visualize the new clusters via the UMAP saved in the Seurat object.

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

The clusters generate by fcoex match possible matches different Seurat clusters. Looking at the HN clusters: M1 matches cluster 1 (likely monocytes), M2 and M4 match clusters 1 and 2 (likely APCs, B + monocytes), M5 matches cluster 2 (likeky B) M7 maches a subset of cluster 0, and as it includes granzymes and granulolysins, this subset of 0 is likely cytotoxic cells (either NK or CD8)

Let’s just take a look at the M2 individually:

  fc <- recluster(fc) 

## Detecting clusters for the following modules:

## [1] "S100A8"

## S100A8

## [1] "HLA-DPB1"

## HLA-DPB1

## [1] "TYMP"

## TYMP

## [1] "HLA-DRB1"

## HLA-DRB1

## [1] "MS4A1"

## MS4A1

## [1] "LYZ"

## LYZ

## [1] "CD7"

## CD7

## [1] "IFI6"

## IFI6

## [1] "HLA-DQA2"

## HLA-DQA2

  pbmc_small <- RunUMAP(pbmc_small, dims = 1:10)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 20:04:33 UMAP embedding parameters a = 0.9922 b = 1.112

## 20:04:33 Read 80 rows and found 10 numeric columns

## 20:04:33 Using Annoy for neighbor search, n_neighbors = 30

## 20:04:33 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 20:04:33 Writing NN index file to temp file D:\biocbuild\bbs-3.13-bioc\tmpdir\RtmpgXMT4D\file217413763470
## 20:04:33 Searching Annoy index using 1 thread, search_k = 3000
## 20:04:33 Annoy recall = 100%
## 20:04:34 Commencing smooth kNN distance calibration using 1 thread
## 20:04:34 7 smooth knn distance failures
## 20:04:35 Initializing from normalized Laplacian + noise
## 20:04:35 Commencing optimization for 500 epochs, with 2410 positive edges
## 20:04:36 Optimization finished

  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)

Detecting anticorrelated genes in the modules

As the dataset used here is a small subset of the original, some cells might be in unexpected clusters.

The modules capture also negative correlations. These can be especially interesting, as they point to complementary cell types. Let’s look if the module genes are enriched in the cluster HP ( header positive, in blue) or HN (header negative, in green).

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

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

##                 p_val avg_log2FC pct.1 pct.2    p_val_adj cluster     gene
## TUBB1    1.642052e-03   6.112278 0.225 0.000 3.776720e-01      HN    TUBB1
## GP9      1.642052e-03   5.401909 0.225 0.000 3.776720e-01      HN      GP9
## NGFRAP1  1.642052e-03   4.828998 0.225 0.000 3.776720e-01      HN  NGFRAP1
## HLA-DRB1 2.464537e-15   6.308812 0.975 0.050 5.668436e-13      HP HLA-DRB1
## HLA-DRA  2.981826e-13   4.489918 0.975 0.400 6.858199e-11      HP  HLA-DRA
## HLA-DRB5 1.046771e-12   4.930441 0.875 0.025 2.407574e-10      HP HLA-DRB5
## HLA-DPA1 1.764119e-12   4.544636 0.900 0.150 4.057473e-10      HP HLA-DPA1
## HLA-DQB1 6.587968e-10   4.585457 0.725 0.025 1.515233e-07      HP HLA-DQB1
## HLA-DMB  1.089737e-07   3.727532 0.600 0.025 2.506394e-05      HP  HLA-DMB
## HLA-DQA1 2.018550e-07   4.452150 0.600 0.050 4.642664e-05      HP HLA-DQA1
## HLA-DMA  1.073164e-06   3.216338 0.575 0.050 2.468277e-04      HP  HLA-DMA
## CFP      2.045245e-06   2.701739 0.575 0.050 4.704064e-04      HP      CFP
## LY86     2.536648e-05   3.931002 0.425 0.025 5.834291e-03      HP     LY86
## LGALS2   4.185935e-05   2.697375 0.475 0.050 9.627650e-03      HP   LGALS2
## NFKBIA   1.147847e-04   1.123156 0.800 0.275 2.640049e-02      HP   NFKBIA
## HVCN1    4.187025e-04   5.321594 0.275 0.000 9.630157e-02      HP    HVCN1
## ASGR1    1.642052e-03   3.106563 0.225 0.000 3.776720e-01      HP    ASGR1
##            module
## TUBB1    HLA-DRB1
## GP9      HLA-DRB1
## NGFRAP1  HLA-DRB1
## HLA-DRB1 HLA-DRB1
## HLA-DRA  HLA-DRB1
## HLA-DRB5 HLA-DRB1
## HLA-DPA1 HLA-DRB1
## HLA-DQB1 HLA-DRB1
## HLA-DMB  HLA-DRB1
## HLA-DQA1 HLA-DRB1
## HLA-DMA  HLA-DRB1
## CFP      HLA-DRB1
## LY86     HLA-DRB1
## LGALS2   HLA-DRB1
## NFKBIA   HLA-DRB1
## HVCN1    HLA-DRB1
## ASGR1    HLA-DRB1

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

## Warning in FindMarkers.default(object = data.use, slot = data.slot, counts =
## counts, : No features pass logfc.threshold threshold; returning empty data.frame

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

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