本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/SKoaR1YQvbgOUx1Iv069MA
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.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会给出警告，表型-细胞关联的结果可能不太可靠。
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.github.io/Scissor/vignettes/Scissor_Tutorial.html
