又一篇单基因的单细胞数据挖掘思路
又一篇单基因的单细胞数据挖掘思路
前面在学习单细胞数据挖掘思路时，曾老板给我们演示了 单细胞转录组中单个基因的挖掘思路，见推文：单个基因在单细胞里面如何分析呢？，当时我将老板的视频学习总结了这篇稿子，最近在《单细胞创作者》群里老板又发了一个关于单细胞单基因挖掘的文献作业。我反复对比了前后两个文章，才发现前面那个帖子其实都没有把曾老板的视频中的最重要的精华部分写出来！！！
单个基因在单细胞里面如何分析呢？ 这篇帖子中 探讨了 关键基因 CCL2 在单细胞中的表达，并以其表达与否将细胞分成了 CCL2+ 细胞和 CCL2-细胞，并展示了这两种细胞中不同的其他细胞类型比例，结果认为 CCL2+ 细胞主要在巨噬细胞中表达（占比高），但是这个结论是有问题的！讨论的精华点就在视频最后几分钟，感兴趣的朋友可以再去看看：
我们今天再来学习一下另一篇文献中的单个基因在单细胞转录组中的挖掘方法~
文献背景
这篇文献于2024年7月30日发表在J Transl Med.杂志上，标题为《Targeting tumor-infiltrating CCR8+ regulatory T cells induces antitumor immunity through functional restoration of CD4+ Tconvs and CD8+ T cells in colorectal cancer》。
文章主要针对 CCR8 这个趋化因子受体，发现其选择性表达于肿瘤浸润的调节性T细胞（Tregs），并描述了结直肠癌（CRC）中浸润的CCR8+ Tregs的详细特征：
- CRC 患者和 CRC 皮下小鼠模型的肿瘤中主要呈现Tregs细胞；
- 肿瘤浸润的 Tregs 上优先表达 CCR8，且与较差的生存相关；
- 大多数 CCR8+ Tregs 表达活化标志物并表现出强烈的抑制功能；
- 使用抗 CCR8 抗体治疗有效地抑制了皮下 CRC 肿瘤的生长，通过有效地恢复CD4+ 常规T细胞（CD4+ Tconvs）和CD8+ T细胞的抗肿瘤免疫力
- 结论：针对 CCR8+ Treg 的靶向治疗在CRC治疗中具有潜在应用价值。
数据集简单介绍
要分析的数据集编号为 GSE108989：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE108989，总共12个样本，下载附件中的文件：GSE108989_CRC.TCell.S11138.count.txt.gz
GSE108989_CRC.TCell.S10805.norm.centered.txt.gz 368.5 Mb (ftp)(http) TXT
GSE108989_CRC.TCell.S11138.TPM.txt.gz 351.6 Mb (ftp)(http) TXT
GSE108989_CRC.TCell.S11138.count.txt.gz 69.7 Mb (ftp)(http) TXT
GSE108989_CRC.bulk.S12.TPM.tab.gz 2.1 Mb (ftp)(http) TAB
GSE108989_CRC.bulk.S12.count.tab.gz 571.0 Kb (ftp)(http) TAB
GSE108989_CRC.bulk.exome.mutation.tab.gz 9.0 Mb (ftp)(http) TAB
数据简单预处理
读取下载好的 GSE108989_CRC.TCell.S11138.count.txt.gz，并创建seurat对象：
###
### Create: Jianming Zeng
### Date: 2023-12-31
### Email: jmzeng1314@163.com
### Blog: http://www.bio-info-trainee.com/
### Forum: http://www.biotrainee.com/thread-1376-1-1.html
### CAFS/SUSTC/Eli Lilly/University of Macau
### Update Log: 2023-12-31 First version
### Update Log: 2024-12-09 by juan zhang (492482942@qq.com)
###
rm(list=ls())
options(stringsAsFactors = F)
library(ggsci)
library(dplyr)
library(future)
library(Seurat)
library(clustree)
library(cowplot)
library(data.table)
library(ggplot2)
library(patchwork)
library(stringr)
library(qs)
library(Matrix)
getwd()
# 创建目录
getwd()
gse <- "GSE108989"
dir.create(gse)
# 方式三：
if(T) {
###### step1: 导入数据 ######
ct <- data.table::fread("GSE108989/GSE108989_CRC.TCell.S11138.count.txt.gz",data.table = F)
ct[1:5, 1:5]
dim(ct)
# 基因去重
kp <- !duplicated(ct$symbol)
table(kp)
ct <- ct[kp, ]
table(is.na(ct$symbol))
ct <- na.omit(ct)
rownames(ct) <- ct[,2]
ct <- ct[,-c(1,2)]
ct[1:5, 1:5]
meta <- read.table("GSE108989/meta.txt", header = T,sep = "\t")
dim(meta)
head(meta)
rownames(meta) <- meta$UniqueCell_ID
table(meta$Patient_ID)
table(meta$majorCluster)
table(meta$sampleType)
# 创建对meta# 创建对象
sce.all <- CreateSeuratObject(counts = ct, meta.data = meta, min.cells = 3)
sce.all
}
# 查看特征
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
sce.all$orig.ident <- sce.all$Patient_ID
table(sce.all$orig.ident)
table(sce.all$majorCluster)
library(qs)
qsave(sce.all, file="GSE108989/sce.all.qs")
然后经过简单的质控与harmony，umap如下：
Treg细胞鉴定
Treg细胞鉴定，随便搜了一个稿子：Treg细胞简介及指标选择，其中的细胞marker如下：
选择res=1.0的结果进行注释：
rm(list=ls())
library(Seurat)
library(ggplot2)
library(SCP) # https://zhanghao-njmu.github.io/SCP/index.html
# https://scillus.netlify.app/vignettes/plotting
library(Scillus) # https://mp.weixin.qq.com/s/Z69GmXORqKczTsMQ68D4Vw
# https://samuel-marsh.github.io/scCustomize/
library(scCustomize)
library(qs)
###### step4: 看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
sce.all.int <- qread("2-harmony/sce.all_int.qs")
sce.all.int
head(sce.all.int@meta.data)
temp <- as.data.frame(table(sce.all.int$majorCluster))
table(Idents(sce.all.int))
table(sce.all.int$sampleType)
table(sce.all.int$RNA_snn_res.0.1)
table(sce.all.int$RNA_snn_res.0.3)
table(sce.all.int$RNA_snn_res.0.5)
table(sce.all.int$RNA_snn_res.0.8)
getwd()
dir.create('3-check-by-1')
select_idet <- "RNA_snn_res.1"
sce.all.int$RNA_snn_res.1
sce.all.int <- SetIdent(sce.all.int, value = select_idet)
table(sce.all.int@active.ident)
head(sce.all.int@meta.data)
p <- DimPlot(sce.all.int, reduction = "umap", group.by = select_idet, label = T) +
ggtitle(select_idet)
p
ggsave(plot=p, filename="3-check-by-1/Dimplot_resolution_1-1.pdf",width = 6, height = 6)
# 可视化
# FOXP3
# CD25 对应基因 IL2RA
# CD127 对应基因 IL7R
p <- DotPlot(sce.all.int, features = c("IL2RA","FOXP3","IL7R"), assay='RNA',group.by = "RNA_snn_res.1",cols = c("grey", "red") )
p
FeaturePlot(sce.all.int, features = "FOXP3")
注释过程：cluster0与8基本上就是Treg细胞了
根据CCR8的表达分组
这一次的分组与上一篇 单个基因在单细胞里面如何分析呢？ 有所不同，这里考虑了这个基因表达的密度曲线分布，并选择了 0.8为分界线，一起来看看合不合理：
# 提取Treg cells
table(Idents(sce.all.int))
sce.all.int <- subset(sce.all.int, ident="CD4+ T Treg cells")
sce.all.int
# CCR8 分组
df <- FetchData(object =sce.all.int, vars = c("umap_1", "umap_2", "CCR8"), layer = "data")
head(df)
range(df$CCR8)
ggplot(data = df, aes(CCR8)) +
geom_density() +
geom_vline(xintercept = 0.8, linetype = "dashed", color = "red")
CCR8 基因确实在0.8处出现了一个拐点：
分组：
library(dplyr)
df$Group <- if_else(df$CCR8>0.8, "CCR8+ Treg", "CCR8- Treg")
head(df)
table(df$Group)
# Negative Positive
# 850 974
sce.all.int <- AddMetaData(sce.all.int, metadata = df)
head(sce.all.int)
那这样就得到了 CCR8+ Treg 和 CCR8- Treg细胞了，作业后面的部分就非常容易啦：
【写作任务2】，简简单单差异分析呢？ 来源于文献：读取GSE108989_CRC.TCell.S11138.count.txt.gz走scanpy或者Seurat的降维聚类分群即可，提取里面的treg然后按照目标基因表达量高低分组后差异分析 Chen et al. Journal of Translational Medicine (2024) 22:709 https://doi.org/10.1186/s12967-024-05518-8
文献中的这种思路 你有没有获得灵感呢？
腾讯云开发者
