Scissor：整合scRNAseq和RNAseq数据识别表型相关细胞
Scissor
什么是 Scissor 算法？
Scissor 算法，Single-Cell Identification of Subpopulations with bulk Sample phenOtype coRrelation。此算法用于整合单细胞转录组数据和 bulk 转录组数据，鉴定和表型相关的细胞亚型，本质在于把 bulk 信息"投影"回单细胞水平，增强单细胞的可解释性。这些表型可以为 disease stage, tumor metastasis, treatment response, and survival outcomes。最终，基于表型相关的细胞亚型进一步解析其基因表达和分子功能。
Paper： Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nature Biotechnology (2021). https://doi.org/10.1038/s41587-021-01091-3
Github：https://github.com/sunduanchen/Scissor
算法原理如下：
图1. Scissor 的分析流程
分析和代码
如下我提供了一个封装好的代码（main.R）：
scissor 算法支持二分类表型和生存表型输入，因此封装函数提供了两类模型，第一类为分类模型，第二类为生存模型。调用 Run_scissor_v1.R 函数后运行如下主函数即可：
# 请运行此脚本即可
# 此脚本为主流程
rm(list = ls())
source("./Code//Run_scissor_v1.R") # 调用封装函数
# 第一种 分类模型，会根据输入数据自动识别
#--------- classification------------------------------------------------------#
sce = readRDS("./Data/input/Result_classification/sce.rds")
bulk_exp = readRDS("./Data/input/Result_classification/bulk_exp.rds")
metadata = readRDS("./Data/input/Result_classification/metadata.rds")
res_dir = "./Data/result/Result_scissor_classification" # 结果文件夹
cutoff = 0.2 #这个参数通常代表在算法中用于过滤的细胞比例或阈值，默认为0.2
alpha = 0.01 # 用于平衡 L1 范数和基于网络的惩罚效应的参数；
name = "LUAD" # 名字
org = "hsa" # 物种 人为hsa 小鼠为mmu
Run_scissor_v1(sce = sce,bulk_exp = bulk_exp,metadata = metadata,
res_dir = res_dir,cutoff = cutoff,alpha = alpha,name = name,
org = org)
# 第二种 预后模型，会根据输入数据自动识别
#--------- prognosis-----------------------------------------------------------#
sce = readRDS("./input/Result_prognosis/sce.rds")
bulk_exp = readRDS("./input/Result_prognosis//bulk_exp.rds")
metadata = readRDS("input/Result_prognosis//metadata.rds")
res_dir = "Result_scissor_prognosis"
cutoff = 0.2
alpha = 0.01
name = "LUAD"
org = "hsa"
Run_scissor_v1(sce = sce,bulk_exp = bulk_exp,metadata = metadata,
res_dir = res_dir,cutoff = cutoff,alpha = alpha,name = name,
org = org)
#------------------------------------------------------------------------------#
需要注意的参数：
sce: Seurat 分析结果对象，需要注释好的结果；
bulk_exp : bulk 转录组表达矩阵可以为 RNAseq 或者基因表达芯片数据；
metadata: 表型信息，分类模型列名为：samples status；预后模型列名为：samples times status；通过这些列名可以对使用的模型进行自动判断；
cutoff: 这个参数通常代表在算法中用于过滤的细胞比例或阈值，默认为0.2；
alpha: 用于平衡 L1 范数和基于网络的惩罚效应的参数；越大选取的细胞越少，越小选取的细胞越多（偏向细胞群）。
封装的函数如下（会在主函数中导入，无需运行），仅供学习。
Run_scissor_v1 = function(sce = sce,bulk_exp = bulk_exp,metadata = metadata,
res_dir = res_dir,cutoff = cutoff,alpha = alpha,
name = name,org = org){
suppressPackageStartupMessages({
if (!require(devtools)){
install.packages("devtools")
}
if (!require(BiocManager)){
install.packages("BiocManager")
}
if (!require(Scissor)) {
devtools::install_github('sunduanchen/Scissor')
}
if (!require(Seurat)) {
install.packages("Seurat")
}
if (!require(tidyverse)) {
install.packages("tidyverse")
}
if (!require(optparse)) {
install.packages("optparse")
}
if (!require(data.table)) {
install.packages("data.table")
}
if (!require(clusterProfiler)) {
BiocManager::install("clusterProfiler")
}
if (!require(org.Hs.eg.db)) {
BiocManager::install("org.Hs.eg.db")
}
if (!require(scEasy)) {
devtools::install_github('Farewellznm/scEasy')
}
library(Scissor)
library(Seurat)
library(tidyverse)
library(clusterProfiler)
library(optparse)
library(data.table)
library(org.Hs.eg.db)
})
if (!dir.exists(res_dir)) {
dir.create(res_dir)
}
result_dir <- res_dir
source("./Code/scissor_v5.R")
#------------------01_prepare data------------------------------------
message("01_prepare data")
# load seurat object
sc_dataset <- sce
# load bulk-data
bulk_dataset <- as.matrix(bulk_exp)
# load meta data
metadata <- metadata
bulk_dataset <- bulk_dataset[,metadata$samples]
#------------------02_Run_scissor------------------------------------
message("02_Run_scissor")
if (identical(metadata$samples,colnames(bulk_dataset))) {
if (colnames(metadata)[2] == "time") {
print("Now run survival model")
phenotype <- metadata[,2:3]
colnames(phenotype) <- c("time", "status")
infos <- scissor_v5(bulk_dataset, sc_dataset, phenotype, alpha = alpha,
family = "cox")
# umap figures
Scissor_select <- rep("none", ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos$Scissor_pos] <- "Scissor+"
Scissor_select[infos$Scissor_neg] <- "Scissor-"
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
sc_dataset$scissor <- factor(sc_dataset$scissor,levels = c("Scissor+","Scissor-","none"))
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c("indianred1","royalblue","grey"), pt.size = 1)
ggsave(paste(result_dir,"umap_scissor.png",sep = "/"),width = 7,height = 7)
ggsave(paste(result_dir,"umap_scissor.pdf",sep = "/"),width = 7,height = 7)
Scissor_select <- rep("none", ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos$Scissor_pos] <- "UnFavorable"
Scissor_select[infos$Scissor_neg] <- "Favorable"
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "prognosis")
sc_dataset$prognosis <- factor(sc_dataset$prognosis,levels = c("UnFavorable","Favorable","none"))
DimPlot(sc_dataset, reduction = 'umap', group.by = 'prognosis', cols = c("indianred1","royalblue","grey"), pt.size = 1)
ggsave(paste(result_dir,"umap_scissor_prognosis.png",sep = "/"))
ggsave(paste(result_dir,"umap_scissor_prognosis.pdf",sep = "/"))
## scissor cells
sce_result <- subset(sc_dataset,scissor %in% c("Scissor+","Scissor-"))
saveRDS(sce_result,paste(result_dir,"sce_result_subset_scissor.rds",sep = "/"))
saveRDS(sc_dataset,paste(result_dir,"sce_result_all.rds",sep = "/"))
### T test DEGs
DefaultAssay(sce_result) <- "RNA"
sce_result$scissor = as.vector(sce_result$scissor)
table(sce_result$scissor)
Idents(sc_dataset) <- "scissor"
All_DEG <- Seurat::FindMarkers(sc_dataset,ident.1 = "Scissor+",ident.2 = "Scissor-",logfc.threshold = 0,min.pct = 0.01)
DEG <- All_DEG[which(All_DEG$p_val_adj <= 0.05 & abs(All_DEG$avg_log2FC) > 1),]
write.csv(DEG,paste(result_dir,"DEG_result.csv",sep = "/"))
write.csv(All_DEG,paste(result_dir,"AllDEGenes_result.csv",sep = "/"))
library(scEasy)
dir_org = getwd()
setwd(res_dir)
scEasy::Run_Gene_Enrichment(geneList = rownames(DEG),org = org)
setwd(dir_org)
} else if (colnames(metadata)[2] == "status") {
print("Now run Classification Model")
phenotype <- metadata$status
names(phenotype) <- metadata$samples
tag <- c(0:(length(unique(metadata$status))-1))
infos <- scissor_v5(bulk_dataset, sc_dataset, phenotype, tag = tag, alpha = alpha,
family = "binomial")
#---plot
# umap figures
Scissor_select <- rep("none", ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos$Scissor_pos] <- "Scissor+"
Scissor_select[infos$Scissor_neg] <- "Scissor-"
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
sc_dataset$scissor <- factor(sc_dataset$scissor,levels = c("Scissor+","Scissor-","none"))
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c("indianred1","royalblue","grey"), pt.size = 1)
ggsave(paste(result_dir,"umap_scissor.png",sep = "/"),width = 7,height = 7)
ggsave(paste(result_dir,"umap_scissor.pdf",sep = "/"),width = 7,height = 7)
Scissor_select <- rep("none", ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos$Scissor_pos] <- paste(name,"+",sep = "")
Scissor_select[infos$Scissor_neg] <- paste(name,"-",sep = "")
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "Condition")
sc_dataset$Condition <- factor(sc_dataset$Condition,levels = c(paste(name,"+",sep = ""),paste(name,"-",sep = ""),"none"))
DimPlot(sc_dataset, reduction = 'umap', group.by = 'Condition', cols = c("indianred1","royalblue","grey"), pt.size = 1)
ggsave(paste(result_dir,"umap_scissor_Condition.png",sep = "/"),width = 7,height = 7)
ggsave(paste(result_dir,"umap_scissor_Condition.pdf",sep = "/"),width = 7,height = 7)
## scissor cells
sce_result <- subset(sc_dataset,scissor %in% c("Scissor+","Scissor-"))
saveRDS(sce_result,paste(result_dir,"sce_result_subset_scissor.rds",sep = "/"))
saveRDS(sc_dataset,paste(result_dir,"sce_result_all.rds",sep = "/"))
### T test DEGs
DefaultAssay(sce_result) <- "RNA"
sce_result$scissor = as.vector(sce_result$scissor)
table(sce_result$scissor)
Idents(sc_dataset) <- "scissor"
All_DEG <- Seurat::FindMarkers(sc_dataset,ident.1 = "Scissor+",ident.2 = "Scissor-",logfc.threshold = 0,min.pct = 0.01)
DEG <- All_DEG[which(All_DEG$p_val_adj <= 0.05 & abs(All_DEG$avg_log2FC) > 1),]
write.csv(DEG,paste(result_dir,"DEG_result.csv",sep = "/"))
write.csv(All_DEG,paste(result_dir,"AllDEGenes_result.csv",sep = "/"))
library(scEasy)
dir_org = getwd()
setwd(res_dir)
scEasy::Run_Gene_Enrichment(geneList = rownames(DEG),org = org)
setwd(dir_org)
} else {
stop("Please check you metadata colnames must be samples time status")
}
} else {
stop("Please check you metadata and bulk data,Ensure consistent sample order!")
}
file.remove("Scissor_inputs.RData")
}
结果解释
如下我展示预后模型的分析结果：其中 scissor+的细胞为 UnFavorable（不利预后相关细胞），scissor-细胞为 Favorable（有利预后相关细胞）。
图2. UMAP 展示 Scissor 阳性和阴性细胞及其对应的表型
为了进一步研究 Scissor+细胞的相关功能，接下来可以基于筛选获得的 Scissor+ vs Scissor- 进行差异分析，对差异基因进行富集分析：结果表明不利预后细胞的差异基因主要富集到 cancer，TNF signaling 等细胞通路。
图3. Scissor+ vs Scissor- 差异基因的GO 和 KEGG 富集分析结果
文献解读
(1) 黑色素瘤免疫治疗相关细胞
为了理解 ICB（免疫检查点阻断）应答背后的机制，在黑色素瘤单细胞 RNA 测序（scRNA-seq）数据集应用了 Scissor 方法，以鉴定与 ICB 应答相关的 T 细胞亚群 [Ref1]。结果如下所示：通过 scissor 算法，在 melanoma T 细胞亚型（cluster num = 6，图 4a）中整合 bulkRNAseq（带有已知道 immunotherapy response）的数据。其中 scissor+细胞为免疫治疗有效（Favorable immunotherapy Response）的相关细胞类型，主要分布在 cluster2 和 cluster3 两类细胞亚型中（图 4b 和 4c）。对 scissor+ vs other 细胞进行差异分析和功能富集分析（图 4d-f），在分子和通路水平证实了其和 ICB 有效反应有关系；基于差异基因构建的免疫治疗信号特征（137 个基因），在独立的 ICB 数据集中评估了效果（图 4g-h）。
值得注意的是在这里黑色素瘤是没有免疫治疗的样本，通过这个分析依旧可以识别和免疫治疗相关的细胞信息。
图4. Scissor identification results on melanoma T cells
(2) 识别耐药相关的肿瘤细胞
为了识别 bladder cancer 中对 cisplatin 耐药相关的肿瘤细胞类型 [Ref2]，使用 scissor 算法整合了 bladder cancer 单细胞图谱（图 5a 和图 5b）中上皮细胞（Epi）亚型（图 5c）和 GDSC（Genomics of Drug Sensitivity in Cancer）中 bulkRNAseq 和 cisplatin 药物敏感和耐药表型数据。识别了和 cisplatin 耐药相关的细胞亚群（主要为 c4 和 c5）（图 5e）。并对其进一步分析发现其 DDRs 评分显著高于其他亚型，并且差异基因富集到了和 Platinum drug resistance 相关的通路（图 5H）。
图5. 识别耐药相关的肿瘤细胞
(3)识别临床预后相关细胞亚群
在 LUAD 和 LUSC 大队列（29 个数据集，556 个样本，128w 个细胞）中不同突变（KRAS，EGFR，STK11）的单细胞图谱中，使用 Scissor 整合 TCGA 带有预后信息的样本，识别和预后相关的细胞亚型 [Ref3]。
图6. 识别耐药相关的肿瘤细胞
QA（问题和答疑）
1. 算法使用相关问题：
QA1：我的单细胞数据仅为单一表型，我可以使用带有表型的 bulk 数据用于分析吗？
完全可以！原始文章已经测试了。
QA2：我如何识别和我表型相关的细胞？Scissor+细胞和 Scissor-细胞的区别？
分类模型：一般来说 scissor+ 对应在 metadata 中对应编码为 1 的表型相关细胞，scissor- 对应在 metadata 中对应编码为 0 的表型相关细胞。
生存模型：在生存分析（survival analysis）中，0 和 1 通常表示的是结局事件（event）的发生状态，也叫"删失指示变量（event indicator）"：1 = 事件已发生（event occurred），例如：死亡、复发、疾病进展等；0 = 事件未发生（censored，删失）。因此 scissor+ 为不利预后相关细胞，scissor-为有利预后相关细胞。
原始论文描述如下：
2.代码使用相关问题：
QA1：代码中提供的 Seuratv4 和 Seuratv5 版本有啥区别？
对于不同版本的 Seurat 包，输入的 Seurat 对象不一样，其他没有区别。
QA2：我如何确保我给的代码一定适合我的数据？
此分析主要用于单细胞转录组数据和 bulk 数据的整合，如果你的数据不符合标准请先预处理。在运行代码的过程中，务必先测试我提供的示例数据后（这样可以确保环境不存在问题），之后将你的数据整理成符合输入格式的形式，再进行进一步分析。
QA3：后续会更新吗？接下来的计划是什么？如何获取示例文件？
会持续更新。目前版本是一个稳定的版本，主要功能为 Scissor 分析，差异分析和富集分析。后续计划优化代码结果输出分类。以及更多的个性化图表！
获取示例文件请关注私信我即可。
参考文献
[Ref1] Sun D, Guan X, Moran AE, Wu LY, Qian DZ, Schedin P, Dai MS, Danilov AV, Alumkal JJ, Adey AC, Spellman PT, Xia Z. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol. 2022;40:527-538.
[Ref2] Li F, Zhang H, Huang Y, Li D, Zheng Z, Xie K, Cao C, Wang Q, Zhao X, Huang Z, Chen S, Chen H, Fan Q, Deng F, Hou L, Deng X, Tan W. Single-cell transcriptome analysis reveals the association between histone lactylation and cisplatin resistance in bladder cancer. Drug Resist Updat. 2024;73:101059.
[Ref3] Salcher S, Sturm G, Horvath L, Untergasser G, Kuempers C, Fotakis G, Panizzolo E, Martowicz A, Trebo M, Pall G, Gamerith G, Sykora M, Augustin F, Schmitz K, Finotello F, Rieder D, Perner S, Sopper S, Wolf D, Pircher A, Trajanoski Z. High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer. Cancer Cell. 2022;40:1503-1520 e1508.
浙公网安备 33010602011771号