In vallotlab/IDclust: Iterative Differential Clustering for single-cell

knitr::opts_chunk$set(
  warning = FALSE, message = FALSE, tidy = TRUE, cache=TRUE
)
options(timeout = max(1000, getOption("timeout")))

In this article, we will show how to run IDclust on a SingleCellExperiment object of a single- cell RNA dataset of the mouse brain from "Joint profiling of histone modifications and transcriptome in single cells from mouse brain,Chenxu Zhu, Yanxiao Zhang, Yang Eric Li, Jacinta Lucero, . Margarita Behrens, Bing Ren, Nature Methods, 2021 Paired-Tag"

Classical analysis of scRNA dataset with ChromSCape

library(IDclust)
library(ChromSCape)

Data

Download, extract & format scRNA of the mouse brain (Zhu et al., 2021) from the GEO portal.

set.seed(47)
# Download dataset
temp = tempfile()
tempdir = tempdir()
download.file("https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE152020&format=file&file=GSE152020%5FPaired%2DTag%5FH3K27ac%5FDNA%5Ffiltered%5Fmatrix%2Etar%2Egz", temp, quiet = TRUE)
untar(temp, exdir = tempdir)

# Download metadata
annot = tempfile()
download.file("http://catlas.org/pairedTag/cellbrowser/Paired-tag/meta.tsv", annot, quiet = TRUE)
metadata = read.table(annot, sep = "\t", header = TRUE)
rownames(metadata) = metadata$Cell_ID
metadata = metadata[which(metadata$Target == "H3K27ac"),]

features = read.table(file.path(tempdir, "04.Paired-Tag_H3K27ac_DNA_filtered_matrix", "bins.tsv"),
                          row.names = NULL, header = F, sep = "\t")[,1, drop = F]
write.table(features, file = file.path(tempdir, "04.Paired-Tag_H3K27ac_DNA_filtered_matrix", "features.tsv"),
                row.names = F, col.names = F, quote = F)

# Create SingleCellExperiment object
out = ChromSCape::read_sparse_matrix(file.path(tempdir, "04.Paired-Tag_H3K27ac_DNA_filtered_matrix"))
out$datamatrix = out$datamatrix[, match(metadata$Cell_ID, gsub(".*_", "", colnames(out$datamatrix)))]
scExp = ChromSCape::create_scExp(out$datamatrix, out$annot_raw)
SummarizedExperiment::colData(scExp) = cbind(SingleCellExperiment::colData(scExp), metadata)

# Subsample cells
scExp = scExp[,sample(ncol(scExp), 5000, replace = F)]

ChromSCape Analysis

We then run a ChromSCape LSI normalization and dimensionality reduction. We can plot the UMAP and color by the cell type.

scExp <- ChromSCape::find_top_features(scExp,n = 100000, keep_others = FALSE, verbose = FALSE)
scExp <- ChromSCape::feature_annotation_scExp(scExp, ref = "mm10")
scExp <- ChromSCape::normalize_scExp(scExp, type = "TFIDF")
scExp <- ChromSCape::reduce_dims_scExp(scExp, dimension_reductions = c("PCA", "UMAP"), n = 10, remove_PC = "Component_1", verbose = F)

ChromSCape::plot_reduced_dim_scExp(scExp, color_by = "Annotation", reduced_dim = "UMAP")

Classical Louvain clustering

We can run a classical Louvain clustering to see the clusters.

scExp <- ChromSCape::find_clusters_louvain_scExp(scExp)

ChromSCape::plot_reduced_dim_scExp(scExp, color_by = "cell_cluster", reduced_dim = "UMAP")

Iterative Differential Clustering

We can now run the Iterative Differential Clustering, that will re-process and re-cluster each cluster iteratively and find subclusters with significant differences between each other.

By default for a SingleCellExperiment object the processing_ChromSCape function is used for re-processing and the differential_ChromSCape function is used to find significant marker genes.

set.seed(47)
output_dir = "~/Tests/IDC_scExp/"
if(!dir.exists(output_dir)) dir.create(output_dir)
scExp = iterative_differential_clustering(
    scExp,
    output_dir = output_dir,
    plotting = FALSE,
    saving = TRUE,
    n_dims = 10,
    dim_red = "PCA",
    vizualization_dim_red = "UMAP",
    processing_function = processing_ChromSCape,
    differential_function = differential_ChromSCape,
    logFC.th = log2(1.5),
    qval.th = 0.01,
    quantile.activation = 0.7,
    min_frac_cell_assigned = 0.1,
    limit = 5,
    limit_by_proportion = NULL,
    starting.resolution = 0.1,
    starting.k = 100,
    resolution = 0.8,
    k = 100,
    verbose = FALSE
)

We can now read in the output 'IDC_summary' object and plot the cluster hierarchies compared to the author clusters. On this plot, each node is a cluster. The colors represent the distribution of author cluster within each cluster. Link between nodes represent a hierarchy in the iteration. The width of the edges is proportional to the number of genes found.

IDC_summary = qs::qread(file.path(output_dir, "IDC_summary.qs"))
plot_cluster_network(scExp, 
                     IDC_summary = IDC_summary,
                     color_by = "Annotation",
                     node_size_factor = 4,
                     legend = FALSE)

A 'IDcluster' column was added to the SingleCellExperiment object, which we can now project the cluster found this way on the UMAP.

ChromSCape::plot_reduced_dim_scExp(scExp, reduced_dim = "UMAP", color_by = "IDcluster", annotate_clusters = F)

We can also plot particular marker genes in the cluster network by changing the 'color_by' parameter to a gene present in the SingleCellExperiment object.

plot_cluster_network(scExp, 
                     IDC_summary = IDC_summary,
                     color_by = "Foxg1", 
                     threshold_to_define_feature_active = 2,
                     legend = FALSE)