【原】跟NBT学Scissor| bulk RNA + scRNA鉴定与目标表型相关的细胞亚群

Nature Biotechnology上发表了题为“Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data”的研究文章。研究团队开发了Scissor算法，可利用大量单细胞数据和表型信息识别与疾病高度相关的细胞亚群。

本文会尝试介绍一下（1）Scissor算法的工作流程和部分文章内容 ，以及（2）使用R实现方法。

一 Scissor算法的工作流程

Scissor的工作流程如图所示

• 图a：首先需要三个数据文件：单细胞表达矩阵、bulk表达矩阵 和 bulk样本的表型文件。

• 图b：其次计算每对细胞和bulk样本的Pearson相关系数构建相关系数矩阵（关键步骤） ，然后 Scissor会通过优化样本表型Y与相关矩阵S的回归模型。

• 图c：由上述优化模型求解的非零系数β用于选择与目标表型相关的细胞亚群，其中Scissor+ 表示所选择的细胞与感兴趣的表型呈正相关，Scissor-为负相关。

  表型信息可以是连续变量、二分向量或临床生存数据，会分别对应不同的回归方法。

• 图d，e ：为控制假阳性后续会有一个可靠性显著性检验以及 差异表达基因分析、功能富集分析和motif分析等下游分析。

以上，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./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./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会给出警告，表型-细胞关联的结果可能不太可靠。

3.4 查看Scissor结果

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"

3.5 Scissor 结果 umap可视化

#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)

可以发现clusters 1 和 3  是与预后差表型 最相关的细胞亚群，还可以做一下差异分析以及富集分析等后续分析。

3.6 优化模型参数

（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./Scissor/vignettes/Scissor_Tutorial.html

参考资料：

[1] https://sunduanchen./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