我们生信技能树的一个学徒在分析一个GEO数据集：GSE141445 的时候，发现并没有线粒体基因表达，但是他看的那篇挖掘此数据的文献《Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RNA sequencing》中却诡异的出现了线粒体质控小提琴图。这是为什么呢？来一起看一看~
数据集编号为 GSE141445：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141445
首先看看 GSE141445 数据背景
这个数据相关的文章于 2021 年 发表在 Nat Cell Biol 上：
Chen S, Zhu G, Yang Y, Wang F et al. Single-cell analysis reveals transcriptomic remodellings in distinct cell types that contribute to human prostate cancer progression. Nat Cell Biol 2021 Jan;23(1):87-98. PMID: 33420488
数据的描述如下：
为了在单细胞水平上理解前列腺癌的异质性，作者从 12 位患者身上收集了 13 份组织样本（12 份原发性样本和 1 份淋巴结转移样本），并进行了单细胞 RNA 测序（scRNA-seq）。经过标准的数据处理和质量过滤，总共获得了 36,424 个细胞：
挖掘此数据的文章
但是到了另一篇挖掘此数据的文章，也就是学徒看的那篇文献，2024年4月份发表在 Aging 上，标题为《Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RNA sequencing》中，我们看到了完全不一样的诡异描述，并且这个图 A 中还出现了跟学徒分析此数据得不到的线粒体基因表达信息：
此文献详细解读版本：《单细胞+bulkRNA分析前列腺癌中巨噬细胞相关基因》
为什么两个分析一样数据集的文章，前后出现了如此大的描述差异呢？你是选择相信原数据文章还是 挖掘它的文献？
上手分析这个数据看看
先去GEO 下载这个数据：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141445。这个页面提供了两个 矩阵，大小还差别巨大：
1、先看看 GSM4203181_data.raw.matrix.txt.gz
先读取数据：
library(data.table)
Sys.time()
raw.data <- fread( 'GSE141445_RAW/GSM4203181_data.raw.matrix.txt.gz', data.table = F)
Sys.time()
dim(raw.data) #
raw.data[1:4,1:4]
这个数据的信息跟它的原文中经过过滤后的细胞数是可以对应上的，总共为36424个细胞：
简单处理并查看是否有线粒体基因，绘制小提琴图
rownames(raw.data)=raw.data[,1]
raw.data=raw.data[,-1]
head(colnames(raw.data))
g=rownames(raw.data)
g[grepl('^MT',g)]
library(AnnoProbe)
ids=annoGene(g,'SYMBOL','human')
head(ids)
sort(table(ids$chr))
ids=ids[ids$biotypes=='protein_coding',]
library(stringr)
metadata <- as.data.frame(str_split(colnames(raw.data),'[-]', simplify = T))
head(metadata)
table(metadata[,2])
metadata=metadata[,c(2,1)]
colnames(metadata)=c('orig.ident' ,'barcode')
metadata$orig.ident=paste0('p',metadata$orig.ident)
table(metadata$orig.ident)
rownames(metadata) = colnames(raw.data)
identical(rownames(metadata),colnames(raw.data))
sce.all <- CreateSeuratObject(counts = raw.data)
sce.all <- AddMetaData(object = sce.all, metadata = metadata)
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident)
sce.all
# 计算线粒体基因比例
mito_genes <- rownames(sce.all)[grep("^MT-", rownames(sce.all),ignore.case = T)]
print(mito_genes) #可能是13个线粒体基因
sce.all <- PercentageFeatureSet(sce.all, features = mito_genes, col.name = "percent_mito")
fivenum(sce.all@meta.data$percent_mito)
# 可视化
p1 <- VlnPlot(sce.all, group.by = "orig.ident", features = c("nFeature_RNA", "nCount_RNA", "percent_mito"),
pt.size = 0, ncol = 3) + NoLegend()
p1
很明显是没有线粒体基因的：
小提琴图：
2、看看 GSM4203181_data.matrix.txt.gz
这个文件要比上面那个 很明显是经过质控后的 数据大了很多，重新走一遍上面的流程代码即可，发现这两个数据没啥区别，细胞数和基因数是一样的：
那挖掘这个数据的文章做了什么呢？心里邪魅一笑...
更多演示，可观看生信技能树的视频号：
写在文末
如果你也想做单细胞转录组数据分析，最好是有自己的计算机资源哦，比如我们的2024的共享服务器交个朋友福利价仍然是800，而且还需要有基本的生物信息学基础，也可以看看我们的 时隔5年，我们的生信技能树VIP学徒继续招生啦。
