广州站的学员跟完了我们的2天线下单细胞数据挖掘授课加上线上3周直播互动练习和答疑后，开始了做自己领域的单细胞转录组数据挖掘，是GSE141445，一个前列腺癌的公共数据集，对应的文章是2021发表的：《Single-cell analysis reveals transcriptomic remodellings in distinct cell types that contribute to human prostate cancer progression》
但是学员发现没办法根据线粒体基因表达量百分比去进行质量控制，但是看到了同样的数据集的数据挖掘 文章：《Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RNA sequencing》，做出来的qc图里面确实是有线粒体含量的。如下所示：
这样就很尴尬，如果2021的GSE141445数据集作者提供的表达量矩阵里面有线粒体，我们给出来的代码没办法识别，那我认打认罚！如果确实没有，那么原作者，以及一系列基于GSE141445的数据挖掘文章都有瑕疵了。
我们复现的代码是 ：
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]
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)
然后，很明显的可以看到确实是没有线粒体基因：
如果是有线粒体基因，它的规则应该是首字母是MT然后加上一个短横线作为分隔符，比如下面的这样的基因信息 ：
> g[grepl('^MT',g)]
[1] "MT-ATP6" "MT-ATP8" "MT-CO1" "MT-CO2" "MT-CO3"
[6] "MT-CYB" "MT-ND1" "MT-ND2" "MT-ND3" "MT-ND4"
[11] "MT-ND4L" "MT-ND5" "MT-ND6" "MT1A" "MT1B"
[16] "MT1E" "MT1F" "MT1G" "MT1H" "MT1HL1"
[21] "MT1M" "MT1X" "MT2A" "MT3" "MTA1"
但是作者给出来的表达量矩阵里面确实是没有这样的基因。
假如作者在表达量矩阵里面有线粒体基因，但是被我们的正则匹配错过了， 理论上我们可以在top50基因看到线粒体基因 ：
在单细胞转录组数据中，线粒体基因和核糖体基因的表达量往往相对较高，这可能由以下几个原因导致：
高表达水平：
线粒体基因和核糖体基因的表达水平本身就很高，因为它们参与了细胞的基本功能。线粒体是细胞的能量工厂，而核糖体负责蛋白质合成，这两个过程对细胞的生存和功能至关重要。
活跃的细胞通常会有更高的线粒体活动和蛋白质合成需求，因此这些基因的表达量会相应增加。
某些细胞类型可能因为其特定的功能和代谢状态而表达更多的线粒体和核糖体基因。
细胞的生理状态，如增殖状态或应激反应，也会影响线粒体和核糖体基因的表达。
测序技术本身可能存在偏差，导致某些类型的基因（如线粒体基因和核糖体基因）更容易被检测到。
线粒体DNA在细胞中的拷贝数通常远高于核DNA，这可能导致线粒体基因的表达量在单细胞数据中显得更高。
在单细胞制备过程中，细胞的损伤或死亡可能导致线粒体释放其内容物，包括线粒体DNA，这可能会增加线粒体基因的测序读段。
至于是否需要根据线粒体和核糖体基因的表达量对单细胞表达量矩阵进行过滤，取决于研究目的和数据分析的具体需求。以下是一些考虑因素：
研究目的：
如果研究的重点是线粒体功能或核糖体生物学，那么保留这些基因可能是必要的。 如果研究目的是探索其他生物学过程，如细胞分化、信号传导或疾病机制，过滤掉这些高表达的基因可能有助于减少噪音，专注于感兴趣的生物学变化。
高表达的线粒体和核糖体基因可能导致数据的分布偏差，影响后续分析的准确性。在这种情况下，过滤可能是必要的。
线粒体基因的表达量可以反映细胞的代谢状态和能量需求。如果这些信息对研究很重要，可能不需要过滤。
不同类型的细胞可能有不同的线粒体和核糖体基因表达模式。在某些情况下，这些模式可能提供有关细胞身份和功能的有用信息。
在进行差异表达分析时，线粒体和核糖体基因的高表达量可能会影响统计测试的结果。过滤这些基因可以减少这种影响。
写在文末
如果你也想做单细胞转录组数据分析，最好是有自己的计算机资源哦，比如我们的2024的共享服务器交个朋友福利价仍然是800，而且还需要有基本的生物信息学基础，也可以看看我们的生物信息学马拉松授课（买一得五） ，你的生物信息学入门课。而且下周六日我们在长沙线下授课哦：千呼万唤，让我们长沙线下约起
