本文记录了：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