单细胞的这种命名方式可取吗？

生信技能树

发布于 2025-03-17 17:00:11

2500

代码可运行

运行总次数：0

代码可运行

今天又来学习一篇单细胞数据文章，标题为《Integration of Single-Cell and Bulk Transcriptomes to Identify a Poor Prognostic Tumor Subgroup to Predict the Prognosis of Patients with Early-stage Lung Adenocarcinoma》，这篇文献于2025年1月21号发表在J Cancer杂志上。

这是曾老板发在群里的一个数据集：

【写作任务1】这个文献不应该是出现一个， SLC45A3+ cells were identified as a group of cells lacking significant features of known cell markers. 降维聚类分群就好奇怪， eight groups: T, NK, B, epithelial, myeloid, stromal, and SLC45A3+ cells 大家试试看，做一下，然后写一个笔记

文献中的注释结果

文章中的单细胞数据分析结果如下图：先不说作者很诡异的 tSNE 和 UMAP都是分群只放一个就行他偏偏要两个都展示一下感觉有凑图的嫌疑吧。

在众多分群中有一个群：SLC45A3+ 命名的群出现感觉有点突兀，而且下面的基因对应的是：SLC5A3。文章中的描述为：

SLC45A3+ cells were identified as a group of cells lacking significant features of known cell markers

那这个群为什么不可能出现在数据中呢？这就来看一下！

这个单细胞对应的数据集为：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117570，为非小细胞性肺癌样本，总共8个样本：

GSM3304007 P1_Tumor
GSM3304008 P1_Normal
GSM3304009 P2_Tumor
GSM3304010 P2_Normal
GSM3304011 P3_Tumor
GSM3304012 P3_Normal
GSM3304013 P4_Tumor
GSM3304014 P4_Normal

数据预处理

去GEO中把 GSE117570_RAW.tar 下载下来，解压，读取：

###
### 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 <- "GSE117570"
dir.create(gse)

# 方式一：标准文件夹
###### step1: 导入数据 ######   
samples <- list.files("GSE117570/",pattern = "processed_data.txt.gz")
samples
scRNAlist <- lapply(samples, function(pro){
#pro <- samples[1]
print(pro)

  ct <- data.table::fread(paste0("GSE117570/",pro),data.table = F)
  ct[1:5, 1:5]
  dim(ct)
  rownames(ct) <- ct[,1]
  ct <- ct[,-1]
  ct[1:5, 1:5]

  sce <- CreateSeuratObject(ct, project=gsub("_processed_data.txt.gz","", pro), min.cells = 3)
return(sce)
})
names(scRNAlist) <-  gsub("_processed_data.txt.gz","", samples)
scRNAlist

# merge
sce.all <- merge(scRNAlist[[1]], y=scRNAlist[-1], add.cell.ids=gsub("_processed_data.txt.gz","", samples))
sce.all <- JoinLayers(sce.all) # seurat v5
sce.all


# 查看特征
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
temp <- str_split(sce.all$orig.ident, pattern = "_", simplify = T)
head(temp)
sce.all$orig.ident <- temp[,1]
sce.all$patient <- temp[,2]
sce.all$group <- temp[,3]
table(sce.all$orig.ident)
table(sce.all$patient)
table(sce.all$group)

library(qs)
qsave(sce.all, file="GSE117570/sce.all.qs")

经过简单质控，并降维聚类分群，选择res 0.5的结果，共得到16个亚群，10627个细胞：

使用已知marker进行注释

简单绘制一张已知基因的气泡图如下：

简单注释结果如下：

################## 注释，读取注释文件
sce.all.int
head(sce.all.int@meta.data)
Idents(sce.all.int) <- "RNA_snn_res.0.5"
Idents(sce.all.int)

temp <- read.table("3-check-by-0.5/anno.txt",sep = "\t")
temp

new.cluster.ids <- temp[,2]
names(new.cluster.ids) <- temp[,1]
new.cluster.ids

sce.all.int <- RenameIdents(sce.all.int, new.cluster.ids)
table(Idents(sce.all.int))

##  新增一列注释
anno <- as.data.frame(sce.all.int@active.ident)
sce.all.int <- AddMetaData(sce.all.int, metadata = anno, col.name = "celltype")
head(sce.all.int@meta.data)
table(sce.all.int$celltype)

# 美化版
p <- CellDimPlot(sce.all.int, group.by = "celltype", reduction = "UMAP", label = T,label.size = 3, label_repel = T, label_insitu = T,
                 label_point_size = 1, label_point_color =NA ,label_segment_color = NA)
p
ggsave(plot=p, filename="3-check-by-0.3/Dimplot_celltype.pdf",width = 8, height = 8)

phe <- sce.all.int@meta.data
head(phe)
saveRDS(phe, file = "3-check-by-0.5/phe.RData")

没有什么大问题，再来看看文章中描述的那一群，当我想看看 SLC45A3 这个基因在哪一群中高表达时：

VlnPlot(sce.all.int, features = c("SLC45A3"))

结果却告诉我，额有点懵：

> VlnPlot(sce.all.int, features = c("SLC45A3"))
Error in `FetchData()`:
! None of the requested variables were found: SLC45A3
Run `rlang::last_trace()` to see where the error occurred.

我还担心是不是出问题了，使用了模糊搜索：

grep("SLC45A3",rownames(sce.all.int),value = T)
# character(0)

grep("SLC4",rownames(sce.all.int),value = T)
# [1] "SLC4A1AP" "SLC4A7"   "SLC41A3"  "SLC44A1"  "SLC43A3"  "SLC43A2"  "SLC44A2"  "SLC41A1"  "SLC40A1"  "SLC44A4"  "SLC4A2"   "SLC48A1"  "SLC46A3"  "SLC44A3"

这就是不可能出现的原因吗？这个数据根本就没有 SLC45A3 这个基因！

再来看看这个SLC5A3基因呢：

VlnPlot(sce.all.int, features = c("SLC5A3"),group.by = select_idet)
DotPlot(sce.all.int, features = c("SLC5A3"), assay='RNA',group.by = select_idet,cols = c("grey", "red") )