Nature Biotechnology上发表了题为“Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data”的研究文章。研究团队开发了Scissor算法,可利用大量单细胞数据和表型信息识别与疾病高度相关的细胞亚群。
本文会尝试介绍一下(1)Scissor算法的工作流程和部分文章内容 ,以及(2)使用R实现方法。
一 Scissor算法的工作流程
Scissor的工作流程如图所示
以上,Scissor就可以利用单细胞表达数据,bulk表达矩阵和表型信息,从单细胞数据中自动识别与给定表型相关度最高的细胞亚群。
二 文章TCGA-LUAD数据集结果
文章使用Scissor在模拟数据集,TCGA肺腺癌数据集,TP53 表型, 黑色素瘤scRNA数据集 , 阿尔茨海默症(AD)数据集的测试中均有很好的表现。
下面简单的介绍下TCGA-LUAD数据集的结果,如图j 所示通过计算得到了361 Scissor+ cell 和534 Scissor - cell,对照图i 发现大部分为cancer cell ,只有零星几个是T cell 和 B cell 。图k给了详细的数值,其中Scissor+ cell 中有98.1%为cancer cell 。
表明在有表型数据和bulk seq数据的情况下,Scissor能够很好地区分肿瘤细胞和正常细胞,从单细胞数据中准确识别与目标表型相关的大部分细胞。
注意Scissor+ 和 Scissor- 是相对的,取决于表型结局中0 和 1的设置。
更多数据集的结果请详细阅读原文(https://www.nature.com/articles/s41587-021-01091-3)。
三 R实现Scissor算法
3.1 安装,加载
下载Scissor包和示例数据
#Scissor安装和加载
devtools::install_github('sunduanchen/Scissor')
library(Scissor)
#加载LUAD 的scRNA-seq count数据
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))
#加载LUAD 的bulk-seq 数据 和 生存数据
load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
load(url(paste0(location, 'TCGA_LUAD_survival.RData')))
3.2 查看数据情况
(1)查看三个数据集的基本情况
dim(sc_dataset)
#[1] 33694 4102
dim(bulk_dataset)
#[1] 56716 471
bulk_dataset[1:4,1:4]
# TCGA-05-4249 TCGA-05-4250 TCGA-05-4382 TCGA-05-4384
#TSPAN6 57.52398 66.4940573 30.905120 35.10434
#TNMD 0.00000 0.1748819 0.000000 0.00000
#DPM1 98.83813 135.2046911 84.362971 89.10553
#SCYL3 16.23339 6.0050820 6.226792 16.93050
head(bulk_survival)
# TCGA_patient_barcode OS_time Status
#1 TCGA-05-4249 1158 0
#2 TCGA-05-4250 121 1
#3 TCGA-05-4382 607 0
#4 TCGA-05-4384 426 0
#5 TCGA-05-4389 1369 0
#6 TCGA-05-4390 1126 0
(2)准备一个seurat数据( 最好包含预处理数据和构造网络数据)
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)
names(sc_dataset)
#[1] "RNA" "RNA_nn" "RNA_snn" "pca" "tsne" "umap"
#umap可视化
DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10)
(3)检查bulk seq数据 和 表型数据
表型文件的第一列为示例id,注意其顺序需要与bulk表达矩阵中的列名顺序相同。此处以生存数据为例,表型文件的第二列是一个二元变量,'1'表示事件(如癌症复发或死亡),'0'表示右截尾 (0,1对应最终的Scissor- Scissor+)。
#检查对应顺序
all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)
#TRUE
phenotype <- bulk_survival[,2:3]
colnames(phenotype) <- c("time", "status")
head(phenotype)
# time status
#1 1158 0
#2 121 1
#3 607 0
#4 426 0
#5 1369 0
#6 1126 0
3.3 运行 Scissor
此处为生存数据,所以family = "cox" ,如果是logistic回归则使用family = "binomial"。
infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05,
family = "cox", Save_file = 'Scissor_LUAD_survival.RData')
可以看到Scissor首先输出相关系数的五分位数的summary,相关系数均为正且值都不接近于0。
如果自己使用的数据集的中值相关性过低(< 0.01),则scissors会给出警告,表型-细胞关联的结果可能不太可靠。
Scissor 鉴定出预后较差的 201 Scissor+ cells 和预后较好的 4个 Scissor- cells ,可以通过infos1数据中查看这205个cell 的具体信息
names(infos1)
#[1] "para" "Coefs" "Scissor_pos" "Scissor_neg"
length(infos1$Scissor_pos)
#[1] 201
infos1$Scissor_pos[1:4]
#[1] "AAAGTAGAGGAGCGAG_19" "AACCATGCATCTCCCA_19" "AACCGCGAGCTGCGAA_20" "AACTCAGTCCGCGGTA_19"
length(infos1$Scissor_neg)
#[1] 4
infos1$Scissor_neg
#[1] "ACGCCAGTCCTCCTAG_20" "ACGGGCTAGTGGCACA_20" "CCGGTAGGTACCCAAT_15" "GACGCGTAGTGGTCCC_20"
#Scissor+ Scissor- 定义
Scissor_select <- rep("Background cells", ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos1$Scissor_pos] <- "Scissor+ cell"
Scissor_select[infos1$Scissor_neg] <- "Scissor- cell"
#metadata 中添加 Scissor信息
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','royalblue','indianred1'),pt.size = 2)
(1)alpha 设置为NULL
infos2 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = NULL, cutoff = 0.03,
family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
(2)alpha 设置梯度
infos3 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = seq(1,10,2)/1000, cutoff = 0.2,
family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
更多参数,模型详见 https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
参考资料:
[1] https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
[2] Sun D, Guan X, Moran AE, et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nat Biotechnol. 2021;10.1038/s41587-021-01091-3