Getting data ready
An example data object
In this package, we include an example data object
sce.example. For more details, please see
help("sce.example"). A quick look at the data object
below.
data(sce.example)
sce.example
#> class: SingleCellExperiment
#> dim: 16497 865
#> metadata(3): cluster_marker_info f_score cluster_order
#> assays(1): logcounts
#> rownames(16497): A1CF A2M ... ZZEF1 ZZZ3
#> rowData names(1): marker_gene
#> colnames(865): A01_1_Nuclei_NeuNP_H200_1025_MTG_layer1_BCH9
#> A01_BCH3_1NeuNP_H200.1030_MTG_Layer_1 ...
#> P11_1_Nuclei_NeuNN_H200_1025_MTG_layer1_BCH9
#> P11_1_Nuclei_NeuNN_H200_1030_MTG_layer1_BCH8
#> colData names(1): cluster_membership
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
The naming convention and data structure are listed below.
Layer 1 and full MTG data
Here, we introduce two datasets to be used in this vignette. Both
datasets are generated by the Allen Institute for Brain Science using
Smart-seq2 protocol for single nucleus RNA sequencing. The first is the
pre-containerized dataset in the data object sce.example,
which is sampled from a single layer (layer 1) of cerebral cortex from
the middle temporal gyrus (MTG) region of human brain (Boldog et al.,
2018). In the layer 1 dataset, there are 16487 genes and 865 cells;
15 cell type clusters are identified and labeled by the authors.
## rename the data object
sce.layer1 <- sce.example
## cell type clusters and cluster sizes
knitr::kable(table(colData(sce.layer1)$cluster_membership), col.names=c("Cluster", "Size"))
| e1_e299_SLC17A7_L5b_Cdh13 |
299 |
| g1_g48_GLI3_Astro_Gja1 |
48 |
| g2_g27_APBB1IP_Micro_Ctss |
27 |
| g3_g18_GPNMB_OPC_Pdgfra |
18 |
| g4_g9_MOG_Oligo_Opalin |
9 |
| i10_i16_TSPAN12_Vip_Mybpc1 |
16 |
| i1_i90_COL5A2_Ndnf_Car4 |
90 |
| i2_i77_LHX6_Sst_Cbln4 |
77 |
| i3_i56_BAGE2_Ndnf_Cxcl14 |
56 |
| i4_i54_MC4R_Ndnf_Cxcl14 |
54 |
| i5_i47_TRPC3_Ndnf_Car4 |
47 |
| i6_i44_GPR149_Vip_Mybpc1 |
44 |
| i7_i31_CLMP_Ndnf_Cxcl14 |
31 |
| i8_i27_SNCG_Vip_Mybpc1 |
27 |
| i9_i22_TAC3_Vip_Mybpc1 |
22 |
The second dataset is sampled from the full laminar depth (layer 1-6)
of the MTG cerebral cortex (Hodge et al.,
2019). The MTG dataset consist of 13945 genes and 15603 cells; 75
cell type clusters are identified and labeled by the authors. We are
going to devote the next section to introduce a function that compiles
the necessary input data files from the MTG dataset into a data object
that can be used in our matching function.
The biological ground truth is that these two
datasets are from neuroanatomical overlapping regions, therefore, we
would expect each layer 1 cell type to be matched with one full MTG cell
type. Good matches should reflect the laminar characteristics of these
cell types, i.e. matches should be found in the upper layer cell types
of the full MTG dataset.
Construct data object using make_data_object()
The following raw data files are here.
## read in pieces of input data - this may take a few minutes
cell_by_gene_expression <- fread("cell_by_gene_expression.csv")
cell_cluster_labels <- fread("cell_cluster_labels.csv")
NSForest_marker_genes <- fread("NSForest_marker_genes.csv")
NSForest_fscores <- fread("NSForest_fscores.csv")
MTG_taxonomy <- fread("MTG_taxonomy.csv")$x #need to be a vector
## unique markers
unique_markers <- unique(NSForest_marker_genes$markerGene)
Use the make_data_object() function to compile the
pieces of input data into a data object.
sce.MTG <- make_data_object(dat = cell_by_gene_expression,
tab = cell_cluster_labels,
markers = unique_markers,
## below are optional input data
cluster_marker_info = NSForest_marker_genes,
f_score = NSForest_fscores,
cluster_order = MTG_taxonomy)
Take a look at the data object for the full MTG.
sce.MTG
#> class: SingleCellExperiment
#> dim: 13945 15603
#> metadata(3): cluster_marker_info f_score cluster_order
#> assays(1): logcounts
#> rownames(13945): A1BG A2M ... ZZEF1 ZZZ3
#> rowData names(1): marker_gene
#> colnames(15603): F1S4_160106_001_B01 F1S4_160106_001_C01 ...
#> F2S4_170405_060_F01 F2S4_170405_060_H01
#> colData names(1): cluster_membership
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
It may cause some confusion on the cell-by-gene or gene-by-cell
orientation of the expression data. To be clear, the raw input data
(.csv or .txt file) for gene expression should have cells in rows and
genes in columns (cell-by-gene). After the input data is read into the
data object (e.g. sce.example), the expression data is stored as genes
in rows and cells in columns (gene-by-cell). We made such arrangement
because conventionally R packages, such as Seurat, deal with
gene-by-cell expression data; and Python packages, such as Scanpy, deal
with cell-by-gene expression data.
FR-Match
Cell type “barcode”
In the FR-Match workflow, we utilize the informative marker genes for
dimensionality reduction. Good marker genes display on-off binary
expression pattern producing, in combination, a unique gene expression
“barcode” for each cell type cluster. We provide the
plot_cluster_by_markers() function for plotting the marker
gene expression patterns for a random set of cells in each cluster,
visualizing the cell type “barcode”. Examples below show that the
“barcodes” for excitatory cell type, glial cell type, and inhibitory
cell type from the layer 1 data are distinctively different using the
NS-Forest marker genes.
plot_cluster_by_markers(sce.example, cluster.name = "e1_e299_SLC17A7_L5b_Cdh13", name.self = "Layer1_")
plot_cluster_by_markers(sce.example, cluster.name = "g1_g48_GLI3_Astro_Gja1", name.self = "Layer1_")
plot_cluster_by_markers(sce.example, cluster.name = "i1_i90_COL5A2_Ndnf_Car4", name.self = "Layer1_")
Cluster sizes
Take a look at the cluster sizes in these two datasets.
plot_clusterSize(sce.layer1, sce.MTG, name.E1 = "Layer1", name.E2 = "MTG")
Note that cluster sizes range widely, which may cause the unbalance
cluster size issue in the FR test.
Run FRmatch()
The FRmatch() function is a wrapper function that takes
in two input data objects in the sce.query = and
sce.ref = arguments. By performing this function, it
carries out our matching workflow in default setting. Key steps are
reported while running. To start, let’s regard sce.layer1
as the query data and sce.MTG as the reference data.
rst.layer1toMTG <- FRmatch(sce.query = sce.layer1, sce.ref = sce.MTG, subsamp.size = 10)
#> * Check query data object.
#> * Check reference data object.
#> * Feature selection: 108 out of 157 reference marker genes are presented in the query experiment.
#> * Filtering small clusters: query and reference clusters with less than 5 cells are not considered.
#> * Comparing 15 query clusters with 75 reference clusters...
#> ** Parallel computing with 12 cores.
#> ** method = cluster2cluster | subsamp.size = 10 | subsamp.iter = 1000
#> * Done parallel and returning results.
A few remarks:
- Note that not all clusters are passed to the matching algorithm. In
the default setting, query and reference clusters with less 10 cells are
filtered out. This is can be tuned using the
filter.size =
argument. It is also optional to filter clusters based on F-beta score
using the filter.fscore = argument.
FRmatch() uses an iterative subsampling scheme to
overcome the unbalanced cluster sizes between the query and reference
clusters when applying the FR test. The subsampling parameters are
subsamp.size =, subsamp.iter =, subsamp.seed =.- Due to the iterative scheme, it may take some computational time. We
provide the
numCores = argument to specific the number of
cores to utilize for parallel computing.
Also, we may want to swap the query and reference data to perform
FR-Match in the other direction. In practice, if we want to match
between two independently conducted experiments (for example, each
experiment may focus on a different specimen region, inducing different
cell type clusters and different set of marker genes), different
directions of matching may lead to quiet different matching results. We
recommend to perform both directions of matching, and conclude with a
consensus matching results from both directions.
rst.MTGtolayer1 <- FRmatch(sce.query = sce.MTG, sce.ref = sce.layer1, subsamp.size = 10)
#> * Check query data object.
#> * Check reference data object.
#> * Feature selection: 30 out of 32 reference marker genes are presented in the query experiment.
#> * Filtering small clusters: query and reference clusters with less than 5 cells are not considered.
#> * Comparing 75 query clusters with 15 reference clusters...
#> ** Parallel computing with 12 cores.
#> ** method = cluster2cluster | subsamp.size = 10 | subsamp.iter = 1000
#> * Done parallel and returning results.
The FRmatch() output is a list of results. The best way
to present these results is to use our customized plotting
functions.
Plot bi-directional matching results
The final matching results from both directions can be combined using
the plot_bi_FRmatch() function. The function takes in the
FRmatch() outputs from both directions, and displays the
two-way match (match found in both directions), one-way match (match
found in either direction), and no match in the following plot.
plot_bi_FRmatch(rst.layer1toMTG, rst.MTGtolayer1, name.E1="Layer1_", name.E2="MTG_")
FR-Match uniquely maps cell types reflecting the overlapping
anatomic regions. Using FR-Match, we mapped most of the layer 1
cell types uniquely to one MTG cell type, i.e. 1-to-1 two-way matches.
These matches precisely reflect the overlapping anatomic regions in
these two independent experiments in that the matched MTG cell types all
have an ‘L1’ layer indicator in their nomenclature. The couple
exceptions of 1-to-many two-way matches may suggest under-partitioning
of the layer 1 cell type; but the multiple matches are next to each
other given the order of the MTG cell types follow the hierarchical
taxonomy in Figure 1c of Hodge et
al. (2019).
Plot one-directional matching results
We also provide the function plot_FRmatch() that takes
in single FRmatch() output with argument
type = "matches" by default, and optionally
type = "padj", to visualize the matching results from one
direction.
## to visualize the one-directional matches
plot_FRmatch(rst.layer1toMTG)
## to visualize the adjusted p-values for each query cluster
plot_FRmatch(rst.layer1toMTG, type = "padj")
A few remarks:
- For the plots above, note that the prefixes of these cluster names
are “query.” and “ref.”, which are coded in our program as default
option. We deliberately leave it that way since: 1) it helps our users
to distinguish the direction of the one-directional matching; and 2) we
sort of regard the one-directional matching as an internal step and
strongly recommend to perform the matching in both directions. These
prefixes can be changed by specifying the
prefix= argument
in the FRmatch() function.
- Users may also obtain the actual values from above plots by
specifying the
return.value = TRUE argument in the
corresponding plotting function.
Look at the matching results in the other direction
## to visualize the one-directional matches
plot_FRmatch(rst.MTGtolayer1)
## to visualize the adjusted p-values for each query cluster
plot_FRmatch(rst.MTGtolayer1, type = "padj")
It is interesting to see that the matching performance depends on the
granularity and/or quality of the reference dataset. It is expected that
the reference dataset is more comprehensive, such as the MTG data in
this example, so that the one directional matching results in this
direction align more with the bi-directional matching results.
Other useful functions
Friedman-Rafsky test
We also implemented our own function FRtest() for Friedman-Rafsky
(FR) test with customized options. FR test is a multivariate
generalization of non-parametric two-sample test. It is a graphical
model based on the concept of minimum
spanning tree (MST). The MST provides a way to visualize
high-dimensional clustered data in a low-dimensional visualization. A
minimal working example of FR test and MST visualization is below.
# simulate some synthetic data from the same distribution
samp1 <- matrix(rnorm(1000), nrow = 50) #a 50-by-20 matrix: 50 dimensional, 20 data points
samp2 <- matrix(rnorm(1000), nrow = 50) #a 50-by-20 matrix: 50 dimensional, 20 data points
# FR test with MST plot
FRtest(samp1, samp2, plot.MST = TRUE, main = "Minimum spanning tree")
#> runs runs.samp1 runs.samp2 stat p.value
#> 20.0000000 7.0000000 13.0000000 -0.3378778 0.3677276
In the above test, the p-value suggests that no difference between
the two simulated samples.
We encourage our users to visually examine their interested cell type
clusters on MST plots.
