R Package Scissor 学习笔记

简介
Scissor R包提出的Scissor算法(function Scissor)，这是一种新颖的单细胞数据分析方法。利用批量（bulk）表型从单细胞测序(scRNA-seq)数据中识别与表型高度相关的细胞亚群。
其优点如下:
首先，通过Scissor识别的表型相关的细胞亚群具有独特的分子特性，其中可能涉及到特定表型的关键标记基因和生物学过程。
其次，Scissor不需要对单细胞数据进行任何无监督聚类，这避免了对细胞聚类数或聚类分辨率的主观决策。
最后，Scissor提供了一个灵活的框架，将各种外部表型整合到批量数据（bulk）中，以指导单细胞数据(scRNA-seq)分析，实现在无假设前提下去识别临床和生物学相关细胞的亚群。

Scissor 使用实例
Scissor的输入包括三类数据：单细胞表达矩阵、bulk表达矩阵和感兴趣的表型。
每个bulk的表型注释可以是连续因变量、二元组指标向量或临床生存数据。
在本教程中，我们使用几个在肺腺癌(LUAD) 癌细胞scRNA-seq上的应用作为示例，展示如何在实际应用中执行Scissor。

实例一
在第一个例子中，我们使用带有Scissor生存信息的LUAD bulk样本来识别LUAD癌细胞中的侵袭性癌细胞亚群。
【第一步，准备scRNA-seq数据】

library(Scissor)
##Prepare the scRNA-seq data
#We load in the LUAD scRNA-seq raw counts data
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'scRNA-seq.RData')))#load完以后会自动生成sc_dataset

sc_dataset文件有1.1G，查看一部分，确认一下是什么数据

sc_dataset[1:2,1:4]#在LUAD scRNA-seq数据中，每一行代表一个基因，每一列代表一个癌细胞

数据展示

dim(sc_dataset)

dim()结果

这表明共有33694个基因和4102个癌细胞。
对于Scissor中使用的scRNA-seq数据，我们使用Seurat包中的函数对这些数据进行预处理，并构建一个细胞-细胞相似网络。（Seurat R包包含预处理数据和构造网络的函数）。
在Scissor包中，我们将Seurat分析管道封装到Seurat_preprocessing()函数中

Seurat_preprocessing(
  counts,#一种matrix的对象，具有非标准化的数据，列细胞，行为特征（如基因）
  project = "Scissor_Single_Cell",
  min.cells = 400,#一个特征至少得在多少个细胞中检测到，若要将被排除掉的特征重新引入，就将cutoff的阈值调小一点
  min.features = 0,#一个细胞至少包含多少个特征
  normalization.method = "LogNormalize",#'LogNormalize','CLR','RC'
  scale.factor = 10000,#为细胞水平标准化设置比例因子
  selection.method = "vst",#'vst','mvp','disp'
  resolution = 0.6,#如果要获得更多（更少）的聚类，使其高于（低于）1.0。
  dims_Neighbors = 1:10,#用作输入的缩减维度
  dims_TSNE = 1:10,#将哪些维度用作t-SNE的输入特征。
  dims_UMAP = 1:10,#将哪些维度用作UMAP的输入特征
  verbose = TRUE#打印输出
)

#preprocess this data and construct a cell-cell similarity network
sc_dataset<-Seurat_preprocessing(sc_dataset, verbose = F)

执行这一步骤时报错提示内存爆掉了

内存爆掉提示

解决方法：

#解决内存爆掉的方法
ls()#看work space中有什么变量。
object.size()#看每个变量占多大内存。
memory.size()#查看现在的work space的内存使用
memory.limit()#查看系统规定的内存使用上限。如果现在的内存上限不够用，可以通过memory.limit(newLimit)更改到一个新的上限

解决结果

Seurat_preprocessing()函数输出的数据类型是Seurat object

class(sc_dataset)

sc_dataset的数据类型

它包含了所需的预处理矩阵和构建的细胞-细胞相似性网络，以及其他有用的降维结果，如PCA、t-SNE和UMAP。

names(sc_dataset)

sc_dataset中包含的内容

用户可以在Seurat_preprocessing()函数中调整预处理参数以实现不同的目标，并将其他相关的Seurat函数应用于sc_dataset。
例如，我们可以使用UMAP、tSINE、PCA降维可视化这4102个单元格

DimPlot(sc_dataset, reduction = 'umap', label = T, label.size = 10)

UMAP

DimPlot(sc_dataset, reduction = 'tsne', label = T, label.size = 10)

tSINE

DimPlot(sc_dataset, reduction = 'pca', label = T, label.size = 10)

PCA

PS：用户还可以使用自己的Seurat分析pipeline(需要标准化的数据和构建的网络)或其他工具预处理的scRNA-seq数据集提供Seurat对象。

【第二步，准备bulk数据和表型信息】
PS：bulk表达矩阵标准化后要处理成matrix类型的变量，不能是data.frame

##Prepare the bulk data and phenotype
#load in the preprocessed LUAD bulk expression matrix and the corresponding survival data downloaded from TCGA
location <- "https://xialab.s3-us-west-2.amazonaws.com/Duanchen/Scissor_data/"
load(url(paste0(location, 'TCGA_LUAD_exp1.RData')))#load后自动生成bulk_dataset数据

bulk_dataset文件有217.5M，查看一部分，确认一下是什么数据

bulk_dataset[1:3,1:4]#行是基因，列代表细胞（这里不再是单个细胞，而是病人的bulk样本）

bulk_dataset数据展示

dim(bulk_dataset)

dim()结果

这表明该数据共有56716个基因和471个样本。用户不需要保留单细胞和大样本之间的共同基因，这可以在Scissor中自动实现。
此外，所有这些样本都有临床结果：

load(url(paste0(location, 'TCGA_LUAD_survival.RData')))#load后自动生成bulk_survival数据
head(bulk_survival)

临床结果

bulk_survival是一个三列的data.frame。
第一列是病人的id，其顺序应该与bulk表达量矩阵中的顺序相同。

all(colnames(bulk_dataset) == bulk_survival$TCGA_patient_barcode)#检验一下病人的id顺序是否与bulk表达量矩阵中的顺序相同

检查结果

第二列是OS-time的临床观测值。
第三列是一个二元变量，'1'表示事件(如癌症复发或死亡)，'0'表示右删失（只知道实际生存时间大于观察到的生存时间）。

提取其中的二、三列得到临床表型列表

phenotype <- bulk_survival[,2:3]
colnames(phenotype) <- c("time", "status")
head(phenotype)

临床表型列表

【第三步，执行Scissor，选择与表型高度相关的细胞】

##Execute Scissor to select the informative cells
#use Scissor to select the phenotype-associated cell subpopulations
infos1 <- Scissor(bulk_dataset, sc_dataset, phenotype, alpha = 0.05, 
                  family = "cox", #因为表型信息是临床生存信息，故用COX回归模型
                  Save_file = 'Scissor_LUAD_survival.RData')

Scissor执行结果
如打印消息中所示，Scissor首先输出单细胞和bulk样本之间的相关性。我们发现所有的相关系数都是正的，而且这些值都不接近于0。在实际应用中，如果所使用的数据集的中值相关性过低(< 0.01)，则Scissors将给出警告信息，这可能会导致不可靠的表型-细胞关联。
总的来说，Scissor识别到201个Scissor+细胞与较差的生存相关，以及4个Scissor- 细胞与良好的生存相关，这些结果保存在infos1中
infos1中保持的结果

我们可以通过在Seurat对象中添加新的注释来可视化Scissor选择的细胞

#visualize the Scissor selected cells by adding a new annotation in the Seurat object
Scissor_select <- rep(0, ncol(sc_dataset))#创建一个列表，用来表示4102个细胞，
names(Scissor_select) <- colnames(sc_dataset)#给列表中每一个数赋予细胞编号
Scissor_select[infos1$Scissor_pos] <- 1#被选为Scissor+的细胞赋值为1
Scissor_select[infos1$Scissor_neg] <- 2#被选为Scissor-的细胞赋值为2
sc_dataset <- AddMetaData(sc_dataset, metadata = Scissor_select, col.name = "scissor")#将表示4102个细胞分类的列表添加到sc_dataset这个Seurat对象中
DimPlot(sc_dataset, reduction = 'umap', group.by = 'scissor', cols = c('grey','indianred1','royalblue'), pt.size = 1.2, order = c(2,1))#可视化

可视化结果

【第四步，调整参数】
在上述实现中，我们将参数α设置为0.05 (alpha = 0.05)。参数α平衡了L1规范和基于网络的惩罚的效果。较大的α倾向于强调L1规范以鼓励稀疏性，使得Scissor选择更少的细胞。在应用程序中，Scissor可以自动将回归输入保存到一个RData (Save_file = ' scissor_luad_survivor .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')

Scissor结果

在上面的实现中，我们设置了一个新的百分比截止(cutoff = 0.03)，当所选细胞的总百分比小于3%时停止α搜索。对应的α = 0.2，其中选择78个Scissor+细胞和5个Scissor-细胞。
为了避免在实际应用中任意选择alpha，建议根据选择的细胞占总细胞的百分比来搜索α。
cutoff的默认值为0.2，即Scissor选择的细胞数量不应超过单细胞数据中总细胞的20%。
此外，用户可以自定义α序列和百分比截断，以实现不同的目标。

【第五步，可靠性显著性检验】
可靠性显著性检验的目的是：如果选择的单细胞和bulk数据不适合于表型-细胞关联，那么其相关信息就会较少，也不能很好地与表型标签相关联。因此，相应的预测性能较差，与随机排列的标签区分不明显。

#Reliability significance test
load('Scissor_LUAD_survival.RData')
numbers <- length(infos1$Scissor_pos) + length(infos1$Scissor_neg)#统计被Scissor选择的细胞数目
result1 <- reliability.test(X, Y, network, alpha = 0.05, family = "cox", cell_num = numbers, n = 10, nfold = 10)

我们使用10次交叉验证(nfold = 10)，并将置换时间设置为10 (n = 10)，以节省一些时间。在实际应用中，可靠性显著性检验对于较大的排列时间可能很耗时(默认n = 100)。

检验结果

我们可以看到，Test打印一个测试统计量和一个测试p值。这个p值小于0.05，表明这些关联是可靠的。
输出result1还包含使用真实标签和排列标签计算的评估度量值:

names(result1)

result1包含内容

在我们的研究中，评价指标是线性回归的均方误差(mean squared error, MSE)，分类的ROC曲线下面积(area under the ROC curve, AUC)， Cox回归的一致性指数(concordindex, c-index)。使用真实标签计算出的平均评价测度作为检验统计量。

【第六步，细胞水平的评估】
最后，我们可以用函数evaluate.cell()为每个Scissor选择的细胞获得一些支持信息。

evaluate_summary <- evaluate.cell('Scissor_LUAD_survival.RData', infos1, FDR = 0.05, bootstrap_n = 100)

evaluate.cell（）打印结果

函数evaluate.cell()包含两种策略来评估Scissor选择的每个细胞。
第一种策略关注被选细胞和所有bulk样本的相关值及其对应的p值。

evaluate_summary[1:5,1:4]

第一种策略

我们可以看到，evaluate.cell()报告了一个细胞与所有bulk样本的平均相关性(列Mean correlation)和正/负相关性的百分比(列correlation > 0和correlation < 0)
第二种策略使用非参数bootstrap来评估每个Scissor选择细胞的系数分布：

evaluate_summary[1:5,5:10]

第二种策略

利用bootstrap重采样方法进行了数值计算。evaluate.cell()输出five-number摘要显示每个Scissor的系数范围选定的细胞

实例二
在第二个例子中，我们使用TCGA-LUAD提供的其他表型特征来指导鉴定相同的4102个细胞中的细胞亚群。在这里，我们选择的特征为TP53突变，在人类恶性肿瘤中发现的最常见的抑癌基因突变之一。
【第一步，准备scRAN-seq数据】
因为和实例一中所用数据相同，故此步骤跳过。
【第二步，准备bulk数据和表型数据】

##Prepare the bulk data and phenotype
#We download the TP53 mutation status (mutant or wild-type) from TCGA-LUAD as the phenotypes of bulk samples:
load(url(paste0(location, 'TCGA_LUAD_exp2.RData')))#load完会自动生成bulk_dataset

bulk_dataset文件有229.8M，查看一部分，确认一下是什么数据

bulk_dataset[1:3,1:4]#行是基因，列代表细胞（这里不是单个细胞，而是病人的bulk样本）

bulk_dataset数据展示

dim(bulk_dataset)

dim()结果

这表明该数据共有56716个基因和498个样本。

#加载表型信息
load(url(paste0(location, 'TCGA_LUAD_TP53_mutation.RData')))#load后会自动生成TP53_mutation的列表

table(TP53_mutation)

table()结果

这表明，498个bulk样本中，TP53突变样本有255个，其余为野生型
生成表型列表

phenotype <- TP53_mutation
head(phenotype)

TP53突变表型列表

【第三步，执行Scissor，选择与表型高度相关的细胞】

##Execute Scissor to select the informative cells
#use Scissor to select the phenotype-associated cell subpopulations
#Given the above binary variable with TP53 mutant = 1 and wild-type = 0, we use family = 'binomial' in Scissor to select the phenotype-associated cell subpopulations
tag <- c('wild-type', 'TP53 mutant')
infos4 <- Scissor(bulk_dataset, sc_dataset, phenotype, tag = tag, alpha = 0.5, 
                  family = "binomial", #因为表型信息是二元的'0','1'信息，故用logistic回归模型
                  Save_file = "Scissor_LUAD_TP53_mutation.RData")

Scissor执行结果

PS：不同的0-1编码可能导致相反的解释
【第四步，调整参数】

infos5 <- Scissor(bulk_dataset, sc_dataset, phenotype, tag = tag, alpha = 0.2, 
                 family = "binomial", Load_file = "Scissor_LUAD_TP53_mutation.RData")

Scissor结果

总共鉴定出了与TP53突变体相关的414个Scissor+细胞和与野生型相关的318个Scissor-细胞。
我们可以使用UMAP技降维将这些选定的细胞可视化

Scissor_select <- rep(0, ncol(sc_dataset))
names(Scissor_select) <- colnames(sc_dataset)
Scissor_select[infos5$Scissor_pos] <- 1
Scissor_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))

可视化结果

【第五步，可靠性显著性检验】

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)

result2结果

【第六步，细胞水平评估】
与实例一相似，此处省略