本文记录了：TCellSI这个工具的使用（R语言），平台是Linux。
TCellSI: A novel method for T cell state assessment and its applications in immune environment prediction.
- 在Bulk数据中计算T细胞状态得分；
- 在单细胞中计算T细胞状态得分；
- 单细胞数据创建伪Bulk
安装R包TCellSI
# install.packages("devtools")
devtools::install_github("GuoBioinfoLab/TCellSI")
加载R包TCellSI
library(TCellSI)
一、Bulk data
1. 导入示例数据
示例数据是，log2转换的RNA-seq数据的基因表达数据框，其格式为TPM。
sample_expression <- TCellSI::exampleSample
head(sample_expression)
2. 计算得分
ResultScores <- TCellSI::TCSS_Calculate(sample_expression)
查看结果
ResultScores
每一行对应一个 T 细胞状态，每一列代表一个样本名称。
3. 计算在每个样本中各T细胞状态的所占比例，并可视化
可视化前需将数据格式转变为长数据，并计算比例
data <- ResultScores
data <- data %>% t() %>% as.data.frame
data <- data.frame(sample = rownames(data), data)
data
第一列是样本名，第2列到第9列是8种T细胞状态
修改列名，去掉列名中的下划线，目的是后续画图时图例没有下划线。
names(data)[names(data) == "Progenitor_exhaustion"] <- "Progenitor exhaustion"
names(data)[names(data) == "Terminal_exhaustion"] <- "Terminal exhaustion"
data
转变为长数据
long_data <- melt(data, id.vars = "sample", variable.name = "subcategory", value.name = "value")
head(long_data)
计算在每个样本中各T细胞状态的所占比例。
下面的代码中将T细胞状态那一列转变为因子，是为了后续画图时图例按这个顺序排序。
long_data <- long_data %>%
group_by(sample) %>%
mutate(proportion = value / sum(value)) %>%
ungroup()
long_data$subcategory <- factor(long_data$subcategory,
levels = c(
"Quiescence",
"Helper",
"Proliferation",
"Regulating",
"Cytotoxicity",
"Senescence",
"Progenitor exhaustion",
"Terminal exhaustion"
)
)
head(long_data)
将最终用于画图的数据格式命名为TCellSI_result，并设置配色（与官方配色一致）
TCellSI_result <- long_data
custom_colors <-c(
"Quiescence"= rgb(246, 100, 99, maxColorValue =255),
"Helper"= rgb(55, 157, 165, maxColorValue = 255),
"Proliferation" = rgb(107, 165, 195, maxColorValue = 255),
"Regulating" = rgb(250, 199, 76, maxColorValue = 255),
"Cytotoxicity" = rgb(149, 168, 172, maxColorValue = 255),
"Senescence" = rgb(174, 153, 126, maxColorValue = 255),
"Progenitor exhaustion" = rgb(48, 147, 67, maxColorValue = 255),
"Terminal exhaustion" = rgb(0, 107, 164, maxColorValue = 255)
)
使用ggplot2进行绘图
ggplot(TCellSI_result, aes(x = sample, y = proportion, fill = subcategory)) +
geom_bar(stat = "identity") +
labs(title = "This is title",
x = "",
y = "Relative TCSS (%)") +
scale_y_continuous(labels = scales::percent) +
scale_fill_manual(values = custom_colors) +
theme_classic()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
guides(fill = guide_legend(title = NULL))
TCellSI_result
如果将此方法应用于感兴趣的其他细胞状态，则应编制参考谱，并为细胞状态准备特定的标记基因集。
然后，使用以下函数计算细胞状态的得分。
如果选择不提供参考表达谱而直接进行计算，可以使用参数ref=FALSE，而无需提供参考参数。
OtherScores <- TCellSI::CSS_Calculate(sample_expression, ref=TRUE, reference = XXX, markers = XXX)
参考表达谱和marker的格式如下
#reference
#The self-constructed reference should contain log2-transformed, TPM-normalized gene expression data from RNA-seq or scRNA-seq.
# cell_state1 cell_state2 cell_state3 ...
# DDX11L1 0.32323232 0.54567463 0.32456323
# WASH7P 0.82670591 1.89565638 1.40492732
# MIR6859-1 0.02172025 0.03816506 0.52313432
# MIR1302-2HG 0.00000000 0.00000000 0.00032302
# MIR1302-2 0.00000000 0.00000000 0.00002132
# ...
#markers: A list of multiple cell states containing specific gene sets
#The number of marker genes per cell state can vary.
#$cell_state1
#[1] "XXX" "XXX" "XXX" ...
#$cell_state2
#[1] "XXX" "XXX" "XXX" ...
#$cell_state3
#[1] "XXX" "XXX" "XXX" ...
二、在scRNA-seq数据中使用TCellSI
1. 导入示例数据
示例数据是Seurat包官方教程的数据
pbmc <- readRDS("pbmc.rds")
pbmc
查看这个数据的细胞类型
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
2. 数据格式处理
sample_scRNA <- as.matrix(pbmc@assays$RNA@counts)
sample_scRNA <- log(sample_scRNA + 1)
sample_scRNA[1:5, 1:5]
注：官网上是seurat_obj@assays$RNA@counts，但好像得加个as.matrix()
3. 计算T细胞状态得分
scRNA_scores <- TCellSI::TCSS_scRNAseqCalculate(sample_scRNA, core= 4) # core: default value is 4
dim(scRNA_scores)
scRNA_scores[1:8, 1:5]
注：如果在本地跑，要把core设置为1。本次学习笔记是在云服务器中运行的。
在本地运行默认core = 4时的报错
scRNA_scores <- scRNA_scores %>%
t() %>%
as.data.frame()
head(scRNA_scores)
将T细胞状态得分添加到单细胞seurat对象的meta.data中
pbmc@meta.data$Quiescence <- scRNA_scores$Quiescence
pbmc@meta.data$Regulating <- scRNA_scores$Regulating
pbmc@meta.data$Proliferation <- scRNA_scores$Proliferation
pbmc@meta.data$Helper <- scRNA_scores$Helper
pbmc@meta.data$Cytotoxicity <- scRNA_scores$Cytotoxicity
pbmc@meta.data$Progenitor_exhaustion <- scRNA_scores$Progenitor_exhaustion
pbmc@meta.data$Terminal_exhaustion <- scRNA_scores$Terminal_exhaustion
pbmc@meta.data$Senescence <- scRNA_scores$Senescence
head(pbmc@meta.data)
4. 与VlnPlot()、FeaturePlot()、RidgePlot()和DotPlot()结合进行可视化
VlnPlot(pbmc, features = "Quiescence")
FeaturePlot(object = pbmc, features = "Quiescence")
RidgePlot(pbmc, features = "Quiescence", ncol = 1)
DotPlot(pbmc, features = c("Quiescence", "Regulating", "proliferation", "Helper", "Cytotoxicity", "Progenitor_exhaustion", "Terminal_exhaustion", "Senescence")) + RotatedAxis()
三、用于单细胞数据分析的伪Bulk创建
准备一个表达数据，它是 log2(TPM+1) 或归一化单细胞数据。
在这些数据中，每一行代表一个基因，每一列代表一个细胞 ID（见以下示例）。
还需一个单细胞注释文件，其中包括表达式文件中的细胞注释列和细胞 ID 列（见下面的示例）。
1. 准备数据，处理成相应输入格式
# expression data
# NP710.20180123 NP711.20180123 NP71.20180123 ...
#A1BG 0.070079488 0.216835131 6.269805313
#NAT2 0.002001509 0.003654851 0.003190016
#ADA 0.085464008 0.088970085 0.057264107
#...
# single-cell annotation file
# UniqueCell_ID annotation
# NTH5.20180123 CD4_C01_CCR7
# NTH64.20180123 CD4_C01_CCR7
# NTR57.20180123 CD4_C01_CCR7
# ...
以pbmc为示例数据
表达谱数据如下
sample_scRNA <- as.matrix(pbmc@assays$RNA@counts)
sample_scRNA <- log(sample_scRNA + 1)
sample_scRNA[1:5, 1:5]
注释数据如下
annot <- data.frame(pbmc@active.ident)
annot <- data.frame(UniqueCell_ID = rownames(annot), annotation = annot$pbmc.active.ident)
head(annot)
查看示例数据相关信息
unique(annot$annotation)
length(unique(annot$annotation))
unique(pbmc@meta.data$orig.ident)
9种细胞类型，来自一个样本
2. 得到伪Bulk样本
pseudo_bulk <- TCellSI::create_pseudo_bulk(
annotation_data = annot,
expression_data = sample_scRNA,
cluster_col = "annotation", # the column names of annotation in single-cell annotation file
cell_id_col = "UniqueCell_ID", # the column names of Cell_ID in single-cell annotation file
n_clusters = 9, # number of cell types annotated
factor = 1, # number of samples for downsampling, default is 5
sampling_rate = 0.6 # percentage of cells downsampled, default is 0.6
)
# see examples of the result
# pseudo_bulk, each column represents a newly pseudobulk samples, each row represents a gene.
# CD4_C01_CCR7_bulk CD4_C01_CCR7_bulk.1 CD4_C01_CCR7_bulk.2
#A1BG 0.495165739 0.67542360 0.737122107
#NAT2 0.006033183 0.00337272 0.007104438
#ADA 0.855647562 1.06058830 0.898952625
查看pseudo_bulk
pseudo_bulk[1:5, 1:5]
这步的运行时间
3. 计算T细胞状态得分
Result <- TCSS_Calculate(pseudo_bulk)
这步的运行时间
查看Result
Result
head(colnames(Result))
