Scissor实操（一）：利用Cox回归模型识别特征细胞
前言
在之前的推文：Scissor:整合bulk+scRNA鉴定细胞功能亚群中，Immugent简单介绍了Scissor的开发流程和运用方向，本期推文Immugent将通过代码实操来展示Scissor的具体分析流程。简单来说，Scissor就是将bulk数据来源的一些特征映射到scRNA数据上，从而找出具有特殊意义的细胞群。
下面我们以大家熟悉的TCGA-LUAD数据作为输入数据进行展示，并且使用Cox回归找出对LUAD患者生存有影响的特殊细胞群；后续分析可以基于Scissor挑选出的Scissor+和Scissor-的分组进行下游功能分析，Immugent对后面这部分不做过多介绍，因为不同的科学问题后续分析流程也是不同的。
下面开始对今天的分析流程进行代码展示：
代码展示
分析第一步，准备单细胞输入数据。
#devtools::install_github('sunduanchen/Scissor')
library(Scissor)
#Prepare the scRNA-seq data
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))
dim(sc_dataset)
#[1] 33694 4102
这个示例数据总共有33694个基因和4102个癌细胞。
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F)
class(sc_dataset)
## [1] "Seurat"
## attr(,"package")
## [1] "Seurat"
DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10)
image.png
第二步就需要准备相应的bulk数据和表型特征信息了。
load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))
load(url(paste0(location, 'TCGA_LUAD_survival.RData')))
dim(bulk_dataset)
## [1] 56716 471
all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)
## [1] TRUE
phenotype <- bulk_survival[,2:3]
colnames(phenotype) <- c("time", "status")
head(phenotype)
是我们再熟悉不过的TCGA的数据，包含了各个患者的生存信息。下面就开始执行Scissor来选择特征细胞了，流程很简单，就一行代码。
infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05,
family = "cox", Save_file = 'Scissor_LUAD_survival.RData')
names(infos1)
## [1] "para" "Coefs" "Scissor_pos" "Scissor_neg"
length(infos1$Scissor_pos)
## [1] 201
length(infos1$Scissor_neg)
## [1] 4
从分析结果看，只有很少的一部分细胞符合过滤标准，这样突出了之前Immugent介绍的它存在的缺点之一。
Scissor_select <- rep(0, ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos1$Scissor_pos] <- 1
Scissor_select[infos1$Scissor_neg] <- 2
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','indianred1','royalblue'), pt.size = 1.2, order = c(2,1))
image.png
但是我们可以通过设置不同的cutoff来筛选出不同数量的细胞。
infos2 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = NULL, cutoff = 0.03,
family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
infos3 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = seq(1,10,2)/1000, cutoff = 0.2,
family = "cox", Load_file = 'Scissor_LUAD_survival.RData')
## [1] "alpha = 0.001"
## [1] "Scissor identified 1666 Scissor+ cells and 1870 Scissor- cells."
## [1] "The percentage of selected cell is: 86.202%"
##
## [1] "alpha = 0.003"
## [1] "Scissor identified 1262 Scissor+ cells and 37 Scissor- cells."
## [1] "The percentage of selected cell is: 31.667%"
##
## [1] "alpha = 0.005"
## [1] "Scissor identified 889 Scissor+ cells and 28 Scissor- cells."
## [1] "The percentage of selected cell is: 22.355%"
##
## [1] "alpha = 0.007"
## [1] "Scissor identified 729 Scissor+ cells and 24 Scissor- cells."
## [1] "The percentage of selected cell is: 18.357%"
## [1] "|**************************************************|"
总结
从上面的分析流程我们可以感受到Scissor的一个小bug，不过作者对这种情况给了一种解决办法，我们可以通过调整相应的参数，对模型进行优化。这个过程是很考验分析者的生物学知识的。总的来说，没有固定的参数，需要根据后续结果进行设置。如果阈值卡的越低，挑出的细胞会更加具有代表性，但是细胞数目少；而如果想要得到更多的细胞，就需要牺牲掉一些准确性。
此外，Scissor只能基于某一种表征将细胞分成阳性和阴性，而无法再对细胞进行细分。所以作者在分析LUAD的数据时，使用的是肿瘤细胞而不是免疫细胞的单细胞暑假。因此，如果对主要受单一细胞驱动的表型使用Scissor进行分析，这是非常具有优势的；而对于多细胞驱动的表型，Scissor可能无法捕捉到全部信息。
好啦，本期分享到这就结束啦，我们下期再会！
