摘要：本笔记展示了 Scissors 整合 Bulk 和单细胞测序分析的详细过程。
# 引言
本 Scissor R 包包含了提出的 Scissor 算法（ Scissor 函数），这是一种新颖的单细胞数据分析方法，它利用整体表型来从单细胞测序数据中识别与表型关联最密切的细胞亚群。Scissor 的优势可以简要概括如下：首先，Scissor 识别的与表型相关的细胞表现出独特的分子特性，涉及给定表型的关键标记基因和生物过程。其次，Scissor 不需要对单细胞数据进行无监督聚类，避免了关于细胞簇数量或聚类分辨率的主观决定。最后，Scissor 提供了一个灵活的框架，整合各种外部表型的整体数据来指导单细胞数据分析，实现临床和生物学上相关的细胞亚群的假设自由识别。
在这个教程中，我们使用几个示例来帮助用户在实际应用中执行 Scissor。我们首先加载所需的包：
library(Scissor) |
Loading required package: Seurat
Loading required package: SeuratObject
Loading required package: sp
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
(status 2 uses the sf package in place of rgdal)
Attaching package: 'SeuratObject'
The following objects are masked from 'package:base':
intersect, saveRDS
Loading Seurat v5 beta version
To maintain compatibility with previous workflows, new Seurat objects will use the previous object structure by default
To use new Seurat v5 assays: Please run: options(Seurat.object.assay.version = 'v5')
Loading required package: Matrix
Warning: package 'Matrix' was built under R version 4.3.2
# Scissor 示例
Scissor 输入的三个数据源是：单细胞表达矩阵、整体表达矩阵和感兴趣的表型。整体样本的表型注释可以是连续的因变量、二元组指示向量或临床生存数据。在这个教程中，我们使用几个应用在肺腺癌（LUAD）单细胞 RNA 测序癌细胞上作为示例，展示如何在实际应用中执行 Scissor。
# 使用 Cox 回归应用 Scissor
在我们的第一个示例中，我们使用具有生存信息的 LUAD 整体样本在 Scissor 中识别 LUAD 癌细胞中侵袭性癌细胞亚群。
# 准备单细胞 RNA 测序数据
我们加载 LUAD 单细胞 RNA 测序原始计数数据：
load('scRNA-seq.RData') |
在这个 LUAD 单细胞 RNA 测序数据中，每行代表一个基因，每列代表一个癌细胞。LUAD 单细胞数据的维度是：
dim(sc_dataset) |
[1] 33694 4102
这表明总共有 33,694 个基因和 4,102 个癌细胞。在 Scissor 中使用的单细胞 RNA 测序数据，最好是包含预处理数据和构建的网络的 Seurat 对象。我们使用 Seurat 包中的函数来预处理这些数据并构建细胞 - 细胞相似性网络。在这个包中，我们将 Seurat 分析流程封装到以下函数中：
sc_dataset <- Seurat_preprocessing(sc_dataset, verbose = F) |
Counts matrix provided is not sparse. Creating V5 assay in Seurat Object.
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
Seurat_preprocessing 函数的输出是一个 Seurat 对象：
class(sc_dataset) |
[1] "Seurat"
attr(,"package")
[1] "SeuratObject"
它包含所需的预处理矩阵和构建的细胞 - 细胞相似性网络，以及其他有用的降维结果，如 PCA、t-SNE 和 UMAP。
names(sc_dataset) |
[1] "RNA" "RNA_nn" "RNA_snn" "pca" "tsne" "umap"
用户可以在这个函数中调整预处理参数以实现不同的目标，并应用其他相关的 Seurat 函数到 sc_dataset 。例如，我们可以使用 UMAP 坐标可视化这 4,102 个细胞：
DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10) |
此外，用户还可以提供自己的 Seurat 分析流程的 Seurat 对象（需要标准化数据和构建的网络）或其他工具预处理的单细胞 RNA 测序数据集。
# 准备整体数据和表型
我们加载预处理的 LUAD 整体表达矩阵和从癌症基因组图谱（TCGA）下载的相应生存数据。
load('TCGA_LUAD_exp1.RData') | |
load('TCGA_LUAD_survival.RData') |
与 LUAD 单细胞 RNA 测序数据一样，行和列分别是基因和细胞。LUAD 整体数据集的维度是：
dim(bulk_dataset) |
[1] 56716 471
这显示该数据共有 56,716 个基因和 471 个样本。用户不需要保留单细胞和整体样本之间的共同基因，这可以在 Scissor 中自动完成。此外，所有这些样本都有临床结果：
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
临床生存表型需要一个两列矩阵，列名为”time” 和”status”。第一列包含样本 ID（观察值），应与整体表达矩阵中的列顺序相同。第二列是一个二元变量，其中”1” 表示事件（例如癌症复发或死亡），而”0” 表示右删失。
all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode) |
[1] 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
# 执行 Scissor 以选择信息丰富的细胞
根据以上输入，我们可以使用 Scissor 选择与表型相关的细胞亚群，该亚群由 Cox 回归模型拟合（ family = 'cox' ）：
sc_dataset[["RNA"]] <- as(object = sc_dataset[["RNA"]], Class = "Assay") |
Warning: Assay RNA changing from Assay5 to Assay
infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05, | |
family = "cox", Save_file = 'Scissor_LUAD_survival.RData') |
[1] "|**************************************************|"
[1] "Performing quality-check for the correlations"
[1] "The five-number summary of correlations:"
0% 25% 50% 75% 100%
0.01319134 0.19658357 0.23945759 0.29808510 0.83727807
[1] "|**************************************************|"
[1] "Perform cox regression on the given clinical outcomes:"
[1] "alpha = 0.05"
[1] "Scissor identified 201 Scissor+ cells and 4 Scissor- cells."
[1] "The percentage of selected cell is: 4.998%"
[1] "|**************************************************|"
如打印消息所示，Scissor 首先输出单细胞和整体样本之间相关性的五数概要。我们发现所有相关性都是正的，且这些值不接近 0。在实际应用中，如果使用的数据集的中位数相关性过低（< 0.01），Scissor 会发出警告消息，这可能导致不可靠的表型到细胞关联。
总的来说，Scissor 识别了 201 个与生存状况较差相关的 Scissor+ 细胞和 4 个与良好生存相关的 Scissor- 细胞，保存在 infos1 中：
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"
[4] "AACTCAGTCCGCGGTA_19"
length(infos1$Scissor_neg) |
[1] 4
infos1$Scissor_neg |
[1] "ACGCCAGTCCTCCTAG_20" "ACGGGCTAGTGGCACA_20" "CCGGTAGGTACCCAAT_15"
[4] "GACGCGTAGTGGTCCC_20"
这些选中的细胞占总共 4,102 个细胞的 5%。我们可以通过在 Seurat 对象 sc_dataset 中添加新的注释来可视化 Scissor 选中的细胞：
names(Scissor_select) <- colnames(sc_dataset)Scissor_select[infos1$Scissor_pos] <- 1Scissor_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))
其中 1 代表 Scissor+ 细胞，2 代表 Scissor- 细胞。
# 调整模型参数
在上述实施中，我们将参数 α 设置为 0.05（ alpha = 0.05 ）。参数 α 平衡了 L1 - 范数和基于网络的惩罚效果。较大的 α 倾向于强调 L1 - 范数以促进稀疏性，使得 Scissor 选择的细胞数量少于使用较小 α 的结果。在应用中，Scissor 可以自动将回归输入保存到 RData（ Save_file = 'Scissor_LUAD_survival.RData' ），这对用户调整不同的 α 很方便。我们可以设置 Load_file = 'Scissor_LUAD_survival.RData' 并使用默认的 α 序列（ alpha = NULL ）来运行 Scissor：
infos2 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = NULL, cutoff = 0.03, | |
family = "cox", Load_file = 'Scissor_LUAD_survival.RData') |
[1] "alpha = 0.005"
[1] "Scissor identified 887 Scissor+ cells and 28 Scissor- cells."
[1] "The percentage of selected cell is: 22.306%"
[1] "alpha = 0.01"
[1] "Scissor identified 573 Scissor+ cells and 20 Scissor- cells."
[1] "The percentage of selected cell is: 14.456%"
[1] "alpha = 0.05"
[1] "Scissor identified 201 Scissor+ cells and 4 Scissor- cells."
[1] "The percentage of selected cell is: 4.998%"
[1] "alpha = 0.1"
[1] "Scissor identified 139 Scissor+ cells and 5 Scissor- cells."
[1] "The percentage of selected cell is: 3.510%"
[1] "alpha = 0.2"
[1] "Scissor identified 78 Scissor+ cells and 5 Scissor- cells."
[1] "The percentage of selected cell is: 2.023%"
[1] "|**************************************************|"
在上述实施中，我们设置了一个新的百分比截止值（ cutoff = 0.03 ），当选中细胞的总百分比小于 3% 时停止 α 搜索。相应的 α 等于 0.2，选中了 78 个 Scissor+ 细胞和 5 个 Scissor - 细胞。为了避免在实际应用中任意选择 alpha，我们建议根据选中细胞占总细胞的百分比搜索 α。 cutoff 的默认值为 0.2，即 Scissor 选中的细胞数量不应超过单细胞数据中总细胞数的 20%。此外，用户可以设置自定义的 α 序列和百分比截止值以实现不同的目标。例如，我们可以运行 Scissor 如下：
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 1663 Scissor+ cells and 1870 Scissor- cells."
[1] "The percentage of selected cell is: 86.129%"
[1] "alpha = 0.003"
[1] "Scissor identified 1264 Scissor+ cells and 37 Scissor- cells."
[1] "The percentage of selected cell is: 31.716%"
[1] "alpha = 0.005"
[1] "Scissor identified 887 Scissor+ cells and 28 Scissor- cells."
[1] "The percentage of selected cell is: 22.306%"
[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
在另一个例子中，我们使用 TCGA-LUAD 提供的其他表型特征来指导同一 4102 个癌细胞中细胞亚群体的识别。这里我们关注 TP53，这是在人类恶性肿瘤中最常发现的 mutant 抑癌基因之一。
# 准备整体数据和表型
我们从 TCGA-LUAD 下载了 TP53 mutant 状态（mutant 或 wildtype）作为整体样本的表型：
load('TCGA_LUAD_exp2.RData') | |
load('TCGA_LUAD_TP53_mutation.RData') |
共有 498 个整体样本；其中 255 个检测出 TP53 mutant，其余为 wildtype：
table(TP53_mutation) |
TP53_mutation
0 1
243 255
# 执行 Scissor 以选出信息性细胞
鉴于上述二元变量中 TP53 mutant = 1 和 wildtype = 0 ，我们使用 family = 'binomial' 在 Scissor 中选择与表型相关的细胞亚群体：
phenotype <- TP53_mutation | |
tag <- c('wildtype', 'TP53 mutant') | |
infos4 <- Scissor(bulk_dataset, sc_dataset, phenotype, tag = tag, alpha = 0.5, | |
family = "binomial", Save_file = "Scissor_LUAD_TP53_mutation.RData") |
[1] "|**************************************************|"
[1] "Performing quality-check for the correlations"
[1] "The five-number summary of correlations:"
0% 25% 50% 75% 100%
0.01308976 0.19502200 0.23725380 0.29589213 0.83830571
[1] "|**************************************************|"
[1] "Current phenotype contains 243 wildtype and 255 TP53 mutant samples."
[1] "Perform logistic regression on the given phenotypes:"
[1] "alpha = 0.5"
[1] "Scissor identified 173 Scissor+ cells and 160 Scissor- cells."
[1] "The percentage of selected cell is: 8.118%"
[1] "|**************************************************|"
注意，参数 tag 用于双重检查表型编码。表型中不同的 0-1 编码可能会导致相反的解释。有关更多详细信息，请参阅下一节 Scissor 结果解读。类似地，我们可以通过设置 Load_file = "Scissor_LUAD_TP53_mutation.RData" 来直接调整其他 α 值：
infos5 <- Scissor(bulk_dataset, sc_dataset, phenotype, tag = tag, alpha = 0.2, | |
family = "binomial", Load_file = "Scissor_LUAD_TP53_mutation.RData") |
[1] "alpha = 0.2"
[1] "Scissor identified 414 Scissor+ cells and 320 Scissor- cells."
[1] "The percentage of selected cell is: 17.894%"
[1] "|**************************************************|"
总的来说，Scissor 识别了 414 个与 TP53 mutant 相关的 Scissor+ 细胞和 318 个与 wildtype 相关的 Scissor- 细胞。我们可以使用 UMAP 技术可视化这些选中的细胞：
names(Scissor_select) <- colnames(sc_dataset)Scissor_select[infos5$Scissor_pos] <- 1Scissor_select[infos5$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))
其中 1 代表 Scissor+ 细胞，2 代表 Scissor- 细胞。
# Scissor 结果解读
Scissor+ 细胞和 Scissor- 细胞都是 Scissor 选中的与特定表型高度相关的细胞。这些细胞与表型之间的关联依赖于所使用的模型，并应在特定的上下文中进行解释。对于线性回归和分类模型，参数 phenotype 的初始值会影响解释结果。例如，在上述应用中使用 TP53 mutant 状态，如果将 TP53 mutant 编码为 1 ，wildtype 为 0 ，则 Scissor+ 细胞将与 TP53 mutant 关联，Scissor- 细胞则与 wildtype 关联。如果在 phenotype 参数中对这两种表型的编码进行反转，那么 Scissor+ 细胞和 Scissor- 细胞的解释也会相应改变。对于 Cox 回归，Scissor+ 细胞与较差的生存率相关联，而 Scissor- 细胞则与较好的生存率相关联。
Scissor 可以将细胞与表型关联起来，这种关联是表型之间的相对概念。也就是说，Scissor 指出一个细胞与某个表型相比更可能与另一个表型关联。考虑到单细胞与整体样本之间可能存在的负相关性，我们可以将一个细胞进一步分类为以下三种类别（见表 1）：如果一个细胞与所有整体样本的平均相关性大于 0，并且正相关的数量多于负相关的数量，这个细胞与相关表型更为相似；如果一个细胞的平均相关性小于 0，并且负相关的数量多于正相关的数量，这个细胞应被解释为与另一表型更不相似；否则，这个细胞与表型的关联是不确定的。在大多数情况下，负相关值很少，识别出的细胞归为” 更相似” 类别。
# 可靠性显著性检验
为了确定推断出的表型与细胞关联是否可靠，我们使用 reliability.test 函数进行可靠性显著性检验。进行可靠性显著性检验的动机是：如果选择的单细胞和整体数据不适合表型与细胞关联，则相关性将信息含量较低，与表型标签关联不紧密。因此，相应的预测性能会较差，与随机排列的标签相比无显著区别。在本教程中，我们以上述应用中识别的关联为例，展示如何运行 reliability.test 。
我们选择了与较差或良好生存率相关联的 205 个 Scissor 选中细胞。Scissor 保存的 Rdata 可用作测试输入：
load('Scissor_LUAD_survival.RData') |
我们使用 10 折交叉验证（ nfold = 10 ）并将排列次数设置为 10（ n = 10 ）以节省时间。在实际应用中，对于较大的排列次数（默认 n = 100 ），可靠性显著性检验可能会非常耗时。
numbers <- length(infos1$Scissor_pos) + length(infos1$Scissor_neg) | |
result1 <- reliability.test(X, Y, network, alpha = 0.05, family = "cox", cell_num = numbers, n = 10, nfold = 10) |
[1] "|**************************************************|"
[1] "Perform cross-validation on X with true label"
Finished!
[1] "|**************************************************|"
[1] "Perform cross-validation on X with permutated label"
Finished!
[1] "Test statistic = 0.628"
[1] "Reliability significance test p = 0.000"
我们可以看到 reliability.test 输出了一个测试统计量和一个 p 值。这个 p 值小于 0.05，表明这些关联是可靠的。输出的 result1 还包含了使用真实标签和随机排列标签计算的评估测量值：
names(result1) |
[1] "statistic" "p" "c_index_test_real"
[4] "c_index_test_back"
在我们的研究中，评估测量值是线性回归的均方误差（MSE），分类的接收者操作特征曲线下面积（AUC），以及 Cox 回归的一致性指数（c-index）。使用真实标签计算的平均评估测量值用作测试统计量。
类似地，我们可以测试与 TP53 mutant 或 wildtype 相关联的 Scissor 选中细胞。
load('Scissor_LUAD_TP53_mutation.RData') | |
numbers <- length(infos5$Scissor_pos) + length(infos5$Scissor_neg) | |
result2 <- reliability.test(X, Y, network, alpha = 0.2, family = "binomial", cell_num = numbers, n = 10, nfold = 10) |
Warning: package 'pROC' was built under R version 4.3.2
Type 'citation("pROC")' for a citation.
Attaching package: 'pROC'
The following objects are masked from 'package:stats':
cov, smooth, var
[1] "|**************************************************|"
[1] "Perform cross-validation on X with true label"
Finished!
[1] "|**************************************************|"
[1] "Perform cross-validation on X with permutated label"
Finished!
[1] "Test statistic = 0.790"
[1] "Reliability significance test p = 0.000"
# 细胞水平评估
最后，我们可以使用 evaluate.cell 函数获取每个 Scissor 选中细胞的一些支持信息。以 LUAD 癌细胞的生存应用为例，我们通过运行以下代码评估 205 个 Scissor 选中细胞：
evaluate_summary <- evaluate.cell('Scissor_LUAD_survival.RData', infos1, FDR = 0.05, bootstrap_n = 100) |
[1] "|**************************************************|"
[1] "Performing correlation check for each selected cell"
Finished!
[1] "|**************************************************|"
[1] "Performing nonparametric bootstrap"
Finished!
evaluate.cell 函数包含两种策略来评估每个 Scissor 选中细胞。第一种策略关注每个选中细胞与所有整体样本的相关性值及其对应的 p 值，包括输出的 data.frame 变量中的前四列：
evaluate_summary[1:5,1:4] |
Mean correlation Correlation > 0 Correlation < 0
AAAGTAGAGGAGCGAG_19 0.3056072 100% 0%
AACCATGCATCTCCCA_19 0.1974518 100% 0%
AACCGCGAGCTGCGAA_20 0.3569129 100% 0%
AACTCAGTCCGCGGTA_19 0.3628618 100% 0%
AACTCCCAGCGCCTTG_19 0.2835872 100% 0%
Significant Correlation
AAAGTAGAGGAGCGAG_19 100%
AACCATGCATCTCCCA_19 100%
AACCGCGAGCTGCGAA_20 100%
AACTCAGTCCGCGGTA_19 100%
AACTCCCAGCGCCTTG_19 100%
我们可以看到 evaluate.cell 函数报告了与所有整体样本的平均相关性（列 Mean correlation ）以及正 / 负相关性的百分比（列 Correlation > 0 和 Correlation < 0 ），这与我们提出的一般解释标准中的条件相符合（见表 1）。在这个应用中，所有 Scissor 选中细胞都与相关表型更相似：
all(evaluate_summary$`Mean correlation` & as.numeric(gsub('%','',evaluate_summary$`Correlation > 0`)) > 50) |
[1] TRUE
evaluate.cell 函数还检查了每个选中细胞的标准相关性 p 值得出的 FDR。它输出了每个选中细胞的显著相关性百分比（FDR < 0.05，列 Significant Correlation ）。
第二种策略使用非参数自助法评估每个 Scissor 选中细胞的系数分布：
evaluate_summary[1:5,5:10] |
Coefficient Beta 0% Beta 25% Beta 50% Beta 75%
AAAGTAGAGGAGCGAG_19 0.009203078 0.0006535843 0.006349622 0.009940265 0.01499857
AACCATGCATCTCCCA_19 0.003779788 0.0001131087 0.003247118 0.006729676 0.01058692
AACCGCGAGCTGCGAA_20 0.005839546 0.0001044831 0.003500107 0.007839735 0.01359622
AACTCAGTCCGCGGTA_19 0.004415840 0.0005875274 0.003877294 0.008254801 0.01476508
AACTCCCAGCGCCTTG_19 0.014565920 0.0003173535 0.010897044 0.018264455 0.02771921
Beta 100%
AAAGTAGAGGAGCGAG_19 0.02747730
AACCATGCATCTCCCA_19 0.02185486
AACCGCGAGCTGCGAA_20 0.02852694
AACTCAGTCCGCGGTA_19 0.03463764
AACTCCCAGCGCCTTG_19 0.04282416
通过利用自助重采样， evaluate.cell 函数输出了每个 Scissor 选中细胞的系数五数概要，包括最小值（列 Beta 0% ）、下四分位数（列 Beta 25% ）、中值（列 Beta 50% ）、上四分位数（列 Beta 75% ）和最大值（列 Beta 100% ）。
