simspec：基于细胞簇谱系相似性整合单细胞数据

在单细胞数据分析中，如何对跨不同实验、条件、批次、时间点或其他技术因素的单细胞测序数据进行多样本整合仍然是一个重大的挑战。需要新的计算方法来整合多个样本矫正批次效应，同时保留不同样本之间内在的生物学信息。在这里，我们提出了一种无监督参考的数据表示方法，即细胞簇谱系相似性(CSS)，其中每个细胞由其与不同样本中独立识别的细胞簇的相似性来进行表示。

CSS的概念是用不同样本中细胞簇的相似性来表示每个细胞。由于技术差异可能以不同的方式影响不同的样本，因此不同样本聚类之间的相似性很可能是不可比较的。因此，我们需要对每个样本中不同细胞簇间的相似性分别进行归一化处理，然后将每个细胞对不同样本的归一化相似性合并起来。

原则上，我们可以使用不同的方法计算细胞和细胞簇之间的相似性。同样，也可以使用不同的归一化策略。在simspec包中，我们基于在给定的基因列表(默认是高度变化的基因)中使用Spearman相关性(默认)或Pearson相关性作为相似性的度量。同时，提供了两种可选择的归一化方法: （1）Z-transformation （2）Kernel probability

默认情况下，simspec包使用Spearman相关性和z-transformation归一化方法进行CSS计算。这是对于具有连续分化轨迹和中间状态的发育相关scRNA-seq数据集的首要选择。另一方面，Pearson相关性结合Kernel probability归一化方法更适合用于细胞类型构成较为离散的scRNA-seq数据集。

参考文献：CSS: cluster similarity spectrum integration of single-cell genomics data.

R包安装

# install.packages(devtools)
devtools::install_github("quadbiolab/simspec")

实例演示

接下来，我们将使用不同细胞系的人脑类器官单细胞数据集进行实例演示CSS方法整合单细胞数据。 

数据集：The 10x scRNA-seq data of 20 two-month-old cerebral organoids with four batches (Kanton et al. 2019) 

下载链接：To retrieve the cerebral organoid data, please refer to the data archive in Mendeley Data.

加载所需R包和示例数据集

library(Seurat)
library(simspec)

# counts <- readMM("cerebral_organoids_Kanton_2019_org2m/matrix.mtx.gz")
# meta <- read.table("cerebral_organoids_Kanton_2019_org2m/cells.tsv.gz", sep="\t")
# features <- read.table("cerebral_organoids_Kanton_2019_org2m/features.tsv.gz", sep="\t", stringsAsFactors = F)
# rownames(counts) <- make.unique(features[,2])
# colnames(counts) <- rownames(meta)
# org2m <- CreateSeuratObject(counts = counts, meta.data = meta)

org2m <- readRDS("cerebral_org2m.rds")
org2m
# An object of class Seurat
# 18147 features across 49153 samples within 1 assay
# Active assay: RNA (18147 features, 0 variable features)

head(org2m@meta.data)
#          orig.ident nCount_RNA nFeature_RNA percent.mito batch
#HC00001 SeuratProject      16752         3287   0.01098442     1
#HC00003 SeuratProject      13533         3399   0.01950787     1
#HC00005 SeuratProject       3098         1558   0.02453988     1
 #              organoid   line            barcode region_RSS
#HC00001 h409B2_60d_org1 h409B2 AAACCTGAGAATGTTG-1         MH
#HC00003 h409B2_60d_org1 h409B2 AAACCTGAGAGCCTAG-1         MH
#HC00005 h409B2_60d_org1 h409B2 AAACCTGAGTAATCCC-1         MH
#                                    celltype_RSS region_pred neuron_vs_NPC
#HC00001                   Non-telencephalon NPCs    nontelen    -0.5854159
#HC00003                   Non-telencephalon NPCs    nontelen    -0.7554498
#HC00005 Mesen/rhombencephalon excitatory neurons    nontelen     0.5646920
#          celltype_pred
#HC00001    nontelen_NPC
#HC00003    nontelen_NPC
#HC00005 nontelen_neuron

数据预处理

首先，我们需要对数据进行预处理：数据标准化、筛选高变异基因、PCA降维

org2m <- NormalizeData(org2m)
#Performing log-normalization
#0%   10   20   30   40   50   60   70   80   90   100%
#[----|----|----|----|----|----|----|----|----|----|
#**************************************************|

org2m <- FindVariableFeatures(org2m, nfeatures = 5000)
#Calculating gene variances
#0%   10   20   30   40   50   60   70   80   90   100%
#[----|----|----|----|----|----|----|----|----|----|
#**************************************************|
#Calculating feature variances of standardized and clipped values
#0%   10   20   30   40   50   60   70   80   90   100%
#[----|----|----|----|----|----|----|----|----|----|
#**************************************************|

org2m <- ScaleData(org2m)
#Centering and scaling data matrix
#  |======================================================================| 100%

org2m <- RunPCA(org2m, npcs = 20)
#PC_ 1
#Positive:  STMN2, RTN1, HMP19, MAPT, MLLT11, SOX4, SYT1, GAP43, DCX, GPM6A
#           RAB3A, TMSB10, RBFOX2, NSG1, VAMP2, CD24, SCG3, CRMP1, SEZ6L2, PCSK1N
#           CXADR, MAP1B, GRIA2, STMN4, CNTNAP2, TTC9B, NEUROD2, BEX2, APLP1, RUNX1T1
#Negative:  VIM, HMGB2, SOX2, GNG5, SMC4, ZFP36L1, TUBA1B, SFRP1, HES1, KIAA0101
#           NUSAP1, PHGDH, TTYH1, TYMS, HSPB1, HMGN2, CKS1B, RPS6, FABP7, GSTP1
#           RCN1, CDK1, TOP2A, GAPDH, CKS2, PTTG1, SOX3, HMGN3, BIRC5, MKI67

UMAP降维可视化

org2m <- RunUMAP(org2m, dims = 1:20)
DimPlot(org2m, group.by = "batch") + DimPlot(org2m, group.by = "celltype_RSS")

可以看到，在未进行批次矫正时，不同样本之间相同的细胞类型分散在不同的区域中，存在较大的批次效应。

计算CSS细胞簇谱系相似性

接下来，我们将使用simspec包中的cluster_sim_spectrum函数计算CSS相似性，并进行批次矫正整合不同样本的数据。

label_tag： is the name of a column in the meta.data table that marks the sample/batch information for integration.

org2m <- cluster_sim_spectrum(object = org2m, label_tag = "organoid", 
                              spectrum_type = "corr_ztransform",corr_method = "spearman")

By default, it returns a new Seurat object with an additional reduction object named "css". It can be then used as the reduction methods for following analysis including creating UMAP embedding and clustering. It is worth to mention that all dimensions in the raw CSS output should be used in any following analysis.

org2m <- RunUMAP(org2m, reduction = "css", dims = 1:ncol(Embeddings(org2m, "css")))
org2m <- FindNeighbors(org2m, reduction = "css", dims = 1:ncol(Embeddings(org2m, "css")))
org2m <- FindClusters(org2m, resolution = 1)
DimPlot(org2m, group.by = "batch") + DimPlot(org2m, group.by = "celltype_RSS") + DimPlot(org2m)

可以看到，使用CSS方法进行批次矫正整合后，不同样本之间相同的细胞类型很好的聚集在一起。

比较不同的数据整合方法

接下来，我们将比较CSS方法和常用的一些单细胞数据整合方法：harmony，liger，scanorama，MNN和Seurat CCA等，并使用lisi评估不同数据整合方法的整合效果。

# 加载所需的R包
library(reticulate)                                                       
scanorama <- import("scanorama")
library(Matrix)
library(Seurat)
library(SeuratWrappers)
library(dplyr)
library(simspec)
library(harmony)
library(liger)

library(RANN)
library(lisi)
library(doParallel)

使用不同的数据整合方法进行批次矫正

# integration
## RSS
ref_brainspan <- readRDS("ext/brainspan_fetal.rds")
# 使用RSS方法基于参考数据集进行整合
seurat <- ref_sim_spectrum(seurat, ref_brainspan, reduction.name = "RSS", reduction.key = "RSS_")
seurat <- RunUMAP(seurat, reduction="RSS", reduction.name="umap_RSS", reduction.key="UMAPRSS_", dims = 1:ncol(Embeddings(seurat, "RSS")))

## CSS
seurat <- cluster_sim_spectrum(seurat, label_tag="organoid", cluster_resolution = 0.6)
seurat <- RunUMAP(seurat, reduction="css", dims = 1:ncol(Embeddings(seurat,"css")), reduction.name="umap_css", reduction.key="UMACSS_")

## MNN
seurat_mnn <- RunFastMNN(object.list = SplitObject(seurat, split.by = "organoid"))
seurat_mnn <- RunUMAP(seurat_mnn, reduction="mnn", dims = 1:30)
seurat[['mnn']] <- CreateDimReducObject(Embeddings(seurat_mnn,"mnn")[colnames(seurat),], key="MNN_")
seurat[['umap_mnn']] <- CreateDimReducObject(Embeddings(seurat_mnn,"umap")[colnames(seurat),], key="UMAPMNN_")

## scanorama
mat <- setNames(lapply(SplitObject(seurat, split.by = "organoid"), function(x) as.matrix(t(x@assays$RNA@data[VariableFeatures(seurat),]))), NULL)
gene_list <- lapply(unique(seurat$organoid), function(x) VariableFeatures(seurat))
mat_integrated <- scanorama$integrate(mat, gene_list)
dr_scanorama <- do.call(rbind, mat_integrated[[1]])
rownames(dr_scanorama) <- do.call(c, lapply(mat, rownames))
dr_scanorama <- dr_scanorama[colnames(seurat),]
seurat[['scanorama']] <- CreateDimReducObject(dr_scanorama, key="SCANORAMA_")
seurat <- RunUMAP(seurat, reduction="scanorama", dims = 1:ncol(Embeddings(seurat,"scanorama")), reduction.name="umap_scanorama", reduction.key="UMAPSCANORAMA_")

## Harmony
seurat <- RunHarmony(seurat, "organoid")
seurat <- RunUMAP(seurat, reduction="harmony", reduction.name="umap_harmony", reduction.key="UMAPHARMONY_", dims = 1:50)

## Seurat integration
seurat_samples <- SplitObject(seurat, split.by="organoid")
seurat_samples <- lapply(seurat_samples, function(obj) FindVariableFeatures(obj, nfeatures = 5000))
seurat_anchors <- FindIntegrationAnchors(object.list = seurat_samples, dims = 1:30, anchor.features = 5000)
seurat_integrated <- IntegrateData(anchorset = seurat_anchors, dims = 1:30)
seurat_integrated <- ScaleData(seurat_integrated) %>% RunPCA(npcs = 50, verbose = F) %>% RunUMAP(dims = 1:20)
seurat[['umap_Seurat']] <- CreateDimReducObject(Embeddings(seurat_integrated, "umap_integrated")[colnames(seurat),], key="UMAPSEURAT_")
seurat[['pca_Seurat']] <- CreateDimReducObject(Embeddings(seurat_integrated, "pca")[colnames(seurat),], key = "PCASEURAT_")

## LIGER
seurat_samples <- SplitObject(seurat, split.by="organoid")
data_samples <- lapply(seurat_samples, function(x) x@assays$RNA@data)
data_samples <- lapply(data_samples, function(x) exp(x)-1)
liger_samples <- createLiger(data_samples)
liger_samples <- normalize(liger_samples)
liger_samples <- selectGenes(liger_samples, var.thresh=0.3, do.plot = T)
liger_samples <- scaleNotCenter(liger_samples)
liger_samples <- optimizeALS(liger_samples, k=20, thresh = 5e-5, nrep = 3)
liger_samples <- runUMAP(liger_samples, use.raw=T)
liger_samples <- quantileAlignSNF(liger_samples, resolution = 0.4, small.clust.thresh = 20)
liger_samples <- runUMAP(liger_samples)
coord <- liger_samples@tsne.coords
coord <- coord[colnames(seurat),]
seurat[['umap_LIGER']] <- CreateDimReducObject(coord, key="UMAPLIGER_")
seurat[['LIGER']] <- CreateDimReducObject(liger_samples@H.norm[colnames(seurat),], key="LIGER_")

细胞类型注释

# cell-level annotation
### dorsal-ventral-nontelen
library(glmnet)
library(doMC)
registerDoMC(10)
seurat <- FindNeighbors(seurat, reduction="pca", dims = 1:20) %>%
  FindClusters(resolution = 0.6)

cl_dorsal <- c(2,18,0,12,3,11)
cl_ventral <- c(4,16,7,13,8)
cl_nontelen <- c(19,10,22,20,15)
idx <- c(sample(which(seurat@active.ident %in% cl_dorsal), 1000),
         sample(which(seurat@active.ident %in% cl_ventral), 1000),
         sample(which(seurat@active.ident %in% cl_nontelen), 1000))

x <- as.matrix(seurat@assays$RNA@data[VariableFeatures(seurat),idx])
y <- factor(setNames(rep(c("dorsal","ventral","nontelen"), c(6,5,5)), c(cl_dorsal,cl_ventral,cl_nontelen))[as.character(seurat@active.ident[idx])])
m <- cv.glmnet(x = t(x), y = y, family = "multinomial", parallel = T)
pred_fates <- predict(m, as.matrix(t(seurat@assays$RNA@data[VariableFeatures(seurat),])), type = "response")[,,1]
lab_pred_fates <- factor(colnames(pred_fates)[apply(pred_fates, 1, which.max)])

### NPC-neuron
DE_neuron_NPC <- read.table("ext/DE_NPC_neurons.tsv", sep="\t", stringsAsFactors=F)
high_neuron <- intersect(rownames(seurat), DE_neuron_NPC$feature[which(DE_neuron_NPC$group == "neuron" & DE_neuron_NPC$padj < 0.01 & DE_neuron_NPC$logFC > log(1.2) & DE_neuron_NPC$pct_in > 50 & DE_neuron_NPC$pct_out < 20 & DE_neuron_NPC$auc > 0.6)])
high_NPC <- intersect(rownames(seurat), DE_neuron_NPC$feature[which(DE_neuron_NPC$group == "NPC" & DE_neuron_NPC$padj < 0.01 & DE_neuron_NPC$logFC > log(1.2) & DE_neuron_NPC$pct_in > 50 & DE_neuron_NPC$pct_out < 20 & DE_neuron_NPC$auc > 0.6)])
seurat$neuron_vs_NPC <- colMeans(seurat@assays$RNA@data[high_neuron,]) - colMeans(seurat@assays$RNA@data[high_NPC,])
lab_NPC_neuron <- factor(ifelse(seurat$neuron_vs_NPC > 0, "neuron", "NPC"))

lab_annot <- paste0(lab_pred_fates, "_", lab_NPC_neuron)
seurat$semi_branch <- lab_pred_fates
seurat$semi_celltype <- lab_annot

使用lisi评估不同整合方法的整合效果

# LISI-based benchmark
registerDoParallel(8)
lisi_batch_org2m <- foreach(dr = c("pca","RSS","CSS","mnn","scanorama","harmony","pca_Seurat","LIGER"), .combine = list, .multicombine = T) %dopar%{
  compute_lisi(Embeddings(seurat, dr), seurat@meta.data, "organoid")
}
lisi_ct_org2m <- foreach(dr = c("pca","RSS","CSS","mnn","scanorama","harmony","pca_Seurat","LIGER"), .combine = list, .multicombine = T) %dopar%{
  compute_lisi(Embeddings(seurat, dr), seurat@meta.data, "celltype_reannot")
}
stopImplicitCluster()

layout(matrix(1:2, nrow=1))
boxplot(do.call(cbind, lisi_batch_org2m), frame=F, las=2, names = c("pca","RSS","CSS","mnn","scanorama","harmony","pca_Seurat","LIGER"), outline=F)
boxplot(length(unique(seurat$celltype_reannot))+1-do.call(cbind, lisi_ct_org2m), frame=F, las=2, names = c("pca","RSS","CSS","mnn","scanorama","harmony","pca_Seurat","LIGER"), outline=F)