经常有粉丝在我们单细胞月更群里面问：单个基因在单细胞里面如何分析呢？
一般来说，单个基因的单细胞数据挖掘有两个比较常见的思路：
首先找到任意单细胞数据集然后降维聚类分群后选定一个目标单细胞亚群的特征基因列表去做预后模型； 另外一个思路是：先做预后模型后定位到少量的基因然后去单细胞水平展现它。
我们今天来看看第二种思路，介绍的文献标题为：《Tumor cell-intrinsic MELK enhanced CCL2-dependent immunosuppression to exacerbate hepatocarcinogenesis and confer resistance of HCC to radiotherapy》，于 2024 年 6 月份发表在 Mol Cancer 杂志（IF = 27.7）上。
这篇文章开篇就以 MELK 基因为中心，通过综合运用生物信息学分析、小鼠异种移植实验、HCC肺转移模型和荧光素酶报告基因实验等多种方法，系统地验证了 MELK 在 HCC 中的预后价值及其在肿瘤发生、进展和转移中的关键作用。
1、MELK在HCC中升高，与HCC患者的不良预后相关
TCGA-LIHC队列（图1A）、ICGC队列（图1B）和GSE14520队列（图1C）、68份HCC样本验证集（图1D-F）的分析：与正常组织相比，MELK在肿瘤组织中的表达显著增加，高MELK表达的HCC患者的预后明显差于低MELK表达的患者。
2、接着是各种实验验证：MELK抑制对肝癌发生、发展和自发性肺转移的影响
3、分别找 MELK的上游调控子miRNA，下游调控基因，最后定位到两个关键基因：CCL2、STAT3
MELK的上游调控子miRNA：Mir-505-3p
MELK的下游调控：STAT3/CCL2
4、关键基因 CCL2 在单细胞数据中的展示
这三幅图分别简单的展示了：(P) 在GSE140228队列中，肝细胞癌（HCC）的免疫细胞景观。(Q) 不同免疫细胞中CCL2的表达谱，CD68+巨噬细胞中CCL2表达的增强。 (R) CCL2表达与肿瘤中免疫细胞浸润水平之间的相关性。
GSE140228 这个数据集是非常有名的张泽民院士团队的数据，前面我们也介绍过它：《为什么执着于复现出来张泽民大佬的单细胞分群呢？》。
了解了思路后，我们来简单的复现一下这个图吧，首先是去GEO数据库下载数据，读取并创建seurat对象：
rm(list=ls())
library(ggsci)
library(dplyr)
library(future)
library(Seurat)
library(clustree)
library(Matrix)
getwd()
###### step1: 导入数据 ######
mtx <- readMM( "data/Droplet/GSE140228_UMI_counts_Droplet.mtx.gz" )
cl <- fread( "data/Droplet/GSE140228_UMI_counts_Droplet_barcodes.tsv.gz", header = F,data.table = F )
head(cl)
rl <- fread( "data/Droplet/GSE140228_UMI_counts_Droplet_genes.tsv.gz" ,header = T,data.table = F )
head(rl)
# 基因去重
kp <- !duplicated(rl$SYMBOL)
table(kp)
dim(mtx)
mtx <- mtx[kp,]
rl <- rl[kp,]
rownames(mtx) <- rl$SYMBOL
colnames(mtx) <- cl$V1 #paste0('c',cl$V2)
mtx[1:4,1:4]
dim(mtx)
meta <- fread("data/Droplet//GSE140228_UMI_counts_Droplet_cellinfo.tsv.gz", header = T,data.table = F )
head(meta)
rownames(meta) <- meta[,1]
meta <- meta[,-1]
head(meta)
identical(rownames(meta),colnames(mtx))
library(Seurat)
sce.all <- CreateSeuratObject(counts=mtx, meta.data=meta, min.cells=5, min.features=300)
sce.all
head(sce.all@meta.data)
table(sce.all$Sample)
table(sce.all$Tissue_sub)
table(sce.all$Donor,sce.all$Sample)
table(sce.all$Platform)
table(sce.all$Histology)
table(sce.all$orig.ident)
sce.all$orig.ident <- sce.all$Sample
然后经过简单的质控，harmony以及降维聚类分群，得到以下tsne图：
p <- FeaturePlot(sce.all.filt, features = c("CD68","CCL2"), reduction = "tsne",order = T)
p
以 CCL2基因表达为标准，将细胞亚群分为两大类：Negative 和 Positive，并绘制这两大类中其他细胞亚群的相对比例变化：
# CCL2分组
df <- FetchData(object =sce.all.filt, vars = c("umap_1", "umap_2", "CCL2"), layer = "data")
df$Group <- if_else(df$CCL2>0, "Positive", "Negative")
head(df)
summary(df$CCL2)
table(df$Group)
# Negative Positive
# 64854 476
sce.all.filt <- AddMetaData(sce.all.filt, metadata = df)
head(sce.all.filt)
p <- DimPlot(sce.all.filt, reduction = "tsne", group.by = "Group",label = T,order = T) +
ggtitle("CCL2 Group")
p
使用数据自带的注释绘制两个组别中的细胞亚群相对比例：
library(ggplot2)
dat <- table(sce.all.filt$Group, sce.all.filt$celltype_global)
dat
dat <- dat / rowSums(dat)
dat <- as.data.frame(dat)
dat
colors <- c("#76D7C4", "#F9FB91", "#D3D3D3", "#E74C3C", "#5D8AA8", "#F39C12", "#F8BBD0","#58D68D", "#B3B6B7")
p <- ggplot(dat ,aes(x = Var1, y =Freq, fill = Var2)) +
geom_bar(stat="identity") +
scale_fill_manual(values=colors) +
theme_bw() +
theme(axis.text = element_text(colour = 'black'),
axis.text.x = element_text(angle = 30, hjust = 1)) + #设置主题背景
labs(x = '',y="") #修改x/y轴的标题
p
结果如下：
更多演示，可以观看生信技能树视频号：
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
