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 modules for single cell data. It’s targeted as a tool for exploratory data analysis, and as a stepping stone for obtaining insights from 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 load an a already preprocessed single cell dataset from 10XGenomics. It was 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")

The mini_pbmc3k object is an object of the class SingleCellExperiment that we will explores in the vignette. For more information on the class, check this Introduction to the SingleCellExperiment Class.

Creating the fcoex object

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. (Note: fcoex does not need names, but only clusters, so you can run it even if you do not have names for your clusters.)

# Get clusters from the pre-processing
target <- colData(mini_pbmc3k)
target <- target$clusters

# Get normalized table from the pre-processing
exprs <- as.data.frame(assay(mini_pbmc3k, 'logcounts'))

# Create fcoex object
fc <- new_fcoex(data.frame(exprs),target)

Once you have set up your fcoex object, the first step is to convert a count matrix into a binarized dataframe. The binarization is needed for the fcoex algorithm and empirically we notice that significant biological signal is preserved with the default parameters.

The standard of fcoex 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

Other discretization methods are avaible too, and this step affects the final results in many ways. If you have the time, you can try different discretization parameters and observe the impact in the data.

After the discretization, fcoex builds a co-expression network and extracts modules. Correlations are calculated via a information-theory method called 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.

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] "Number of features features =  1700"
## [1] "Number of prospective features =  199"
## [1] "Number of final features =  28"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

The last step created an adjacency matrix and extracted some co-expression modules. We can present that information in a visual way by using the get_nets function, that generates network plots for each module. The plots are created in the ggplot2 format and stored inside the fcoex object.

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 using the show_net function.

fc <- get_nets(fc)

# Taking a look at the first two networks: 
network_plots <- show_net(fc)

network_plots[["CD79A"]]

network_plots[["HLA-DRB1"]]

Depending on your graphical device, the network might not appear in full. Saving the plots to a vector format, such as pdf, fixes some graphical device issues.

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

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

Additionaly, you can use the ggsave funcion of the ggplot2 package.

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, but you can use any gene set of interest. Be sure, though, that the labels of the genes in the gene sets match the ones in the dataset.

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.

# You'll need CEMiTool, if you do not have it installed, just run BiocManager::install("CEMiTool")
gmt_filename <- system.file("extdata", "pathways.gmt", package = "CEMiTool")

if (gmt_filename == "")
  {
      print("You likely need to install CEMiTool")
} else {
      gmt_in <- pathwayPCA::read_gmt(gmt_filename,  description = TRUE)

}
fc <- mod_ora(fc, gmt_in)

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

## --> Expected input gene ID: RPL8,XRCC5,ABCC1,VCAN,B3GAT1,ABCA7

## --> return NULL...

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

## --> Expected input gene ID: RPS21,EIF3H,PRPS2,RPS17,ABCC9,RPL36A

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

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 pre-processing.

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"

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.

Plotting clusters with schex

After obtaining the new labels, we can visualize them them using UMAP. Let’s see the population represented in the modules CD79A and HLA-DRB1.

For different datasets, different module headers will appear. It is up to the researcher to select which of those are interesting in their research settings. The modules are nevertheless ordered based on correlation with the labels, so the first modules tend to be more interesting.

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.

We will use the package schex for visualizing the data, but you can use your package of preference.

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

identities_based_on_the__CD79A_module <- 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

## Registered S3 method overwritten by 'spatstat.core':
##   method          from
##   formula.glmmPQL MASS

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