线粒体基因缺失之谜

线粒体基因缺失之谜

我们生信技能树的一个学徒在分析一个GEO数据集：GSE141445 的时候，发现并没有线粒体基因表达，但是他看的那篇挖掘此数据的文献《Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RA sequencing》中却诡异的出现了线粒体质控小提琴图。这是为什么呢？来一起看一看~

数据集编号为 GSE141445：.cgi?acc=GSE141445

这个数据相关的文章于 2021 年 发表在 at 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. at Cell Biol 2021 Jan;2(1):87-98. PMID: 420488

数据的描述如下：

为了在单细胞水平上理解前列腺癌的异质性，作者从 12 位患者身上收集了 1 份组织样本（12 份原发性样本和 1 份淋巴结转移样本），并进行了单细胞 RA 测序（scRA-seq）。经过标准的数据处理和质量过滤，总共获得了 6,424 个细胞：

tSE view of 6,424 single cells, color coded by sample (a), cluster (b) and broad lineage (c).

但是到了另一篇挖掘此数据的文章，也就是学徒看的那篇文献，2024年4月份发表在 Aging 上，标题为《Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RA sequencing》中，我们看到了完全不一样的诡异描述，并且这个图 A 中还出现了跟学徒分析此数据得不到的线粒体基因表达信息：

此文献详细解读版本：《单细胞+bulkRA分析前列腺癌中巨噬细胞相关基因》

先去GEO 下载这个数据：.cgi?acc=GSE141445。这个页面提供了两个 矩阵，大小还差别巨大：

1、先看看 GSM420181_data.gz

先读取数据：

library()
()
raw.data <- fread( 'GSE141445_RAW/GSM420181_data.gz',  = F)
()
dim(raw.data) # 
raw.data[1:4,1:4]

这个数据的信息跟它的原文中经过过滤后的细胞数是可以对应上的，总共为6424个细胞：

简单处理并查看是否有线粒体基因，绘制小提琴图

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), = T)] 
print(mito_genes) #可能是1个线粒体基因
sce.all <- PercentageFeatureSet(sce.all, features = mito_genes,  = "percent_mito")
fivenum(sce.all@meta.data$percent_mito)

# 可视化
p1 <- VlnPlot(sce.all, group.by = "orig.ident", features = c("nFeature_RA", "nCount_RA", "percent_mito"), 
              pt.size = 0, ncol = ) + oLegend()
p1 

很明显是没有线粒体基因的：

小提琴图：

2、看看 GSM420181_gz

这个文件要比上面那个 很明显是经过质控后的 数据大了很多，重新走一遍上面的流程代码即可，发现这两个数据没啥区别，细胞数和基因数是一样的：