单细胞层面所谓的无效或者干扰基因去除与否有什么影响呢？

我们生信技能树的一个学徒在看单细胞文献的时候，看到一篇文章在做单细胞数据预处理的时候，过滤了1514个基因！包括50个线粒体基因（通常我们做的时候都是^MT-开头，是13个基因对吧），178个热休克蛋白基因，1253个核糖体基因，以及33个解离相关基因（这个听都没听过啊）。相关的单细胞文献为《UPP1 promotes lung adenocarcinoma progression through the induction of an immunosuppressive microenvironment》，于2024年2月份发表在 Nature Communications杂志上。

下面我们来看看如何获取这1514个基因以及除去与否对结果的影响！

首先，问一下人工智能大模型

问kimi（https://kimi.moonshot.cn/）：人类基因名字里面的：mitochondria(50 genes), heat-shockprotein (178 genes), ribosome (1253genes), and dissociation (33 genes) 这些基因该如何拿到？

要获取人类基因组中与“线粒体（mitochondria）”“热休克蛋白（heat-shock protein）”“核糖体（ribosome）”和“解离（dissociation）”相关的基因信息，可以通过以下途径： 1. 线粒体相关基因（mitochondria genes）

• 数据库查询：可以访问NCBI数据库（National Center for Biotechnology Information），在Gene数据库中搜索“mitochondria genes human”，获取与线粒体相关的基因列表。
• 专业资源：Human MitoCarta 3.0是一个专门的线粒体基因资源库，提供了详细的线粒体相关基因信息，包括基因名称、功能、亚细胞定位等。

2. 热休克蛋白相关基因（heat-shock protein genes）

• 数据库查询：同样可以在NCBI的Gene数据库中搜索“heat-shock protein genes human”，获取相关基因列表。
• 专业文献：根据文献资料，热休克蛋白基因分为多个家族，如HSP33、HSP60、HSP70/HSP110、HSP90、HSP100等。可以参考这些家族名称在数据库中进一步搜索。

3. 核糖体相关基因（ribosome genes）

• 数据库查询：在NCBI的Gene数据库中搜索“ribosome genes human”，可以找到核糖体相关基因。
• 专业文献：文献中提到人类核糖体蛋白基因（RPs）的结构和功能，以及它们在不同物种中的保守性。可以参考这些信息进行更详细的搜索。

4. 解离相关基因（dissociation genes）

• 数据库查询：在NCBI的Gene数据库中搜索“dissociation genes human”，可能会找到与细胞解离或相关功能相关的基因。
• 专业文献：在一些研究中，如单细胞RNA测序（scRNAseq）中提到的细胞解离过程，可能会涉及一些特定基因的表达。可以参考这些研究中的基因列表。

其他建议

• 使用在线工具：一些在线工具如GeneCards（GeneCards）可以提供更全面的基因信息，包括基因功能、表达组织、相关疾病等。
• 参考文献：查阅相关领域的综述文章和研究论文，这些文献中可能会列出与特定功能相关的基因列表。

通过以上方法，可以较为全面地获取这些基因的相关信息。

人工智能大模型这里也没有很智能的给出这些相关基因列表，我们还是多看文献，从文献里面收集。

这些基因其实在另一篇张泽民团队 2022 年发表在 Nature 杂志上的文献附表 table1 中：《Liver tumour immune microenvironment subtypes and neutrophil heterogeneity》：

将这些基因下载并保存在一个txt文档中：

测试数据

我们随便找一个单细胞数据看看去除前与后的情况，这里我们使用 GSE163558 数据：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE163558，这个数据的预处理以及基础分析见我们之前的帖子：《单细胞转录组降维聚类分群过滤基因和过滤细胞的区别》。

1、过滤1514个基因之前的结果

简单看一下数据的降为聚类umap图：

rm(list=ls())
library(Seurat)
library(tidyverse)
library(harmony)

# 加载之前的GSE163558单细胞数据分析结果
data.filt = readRDS('./inputs/sce.all_int.rds')  
data.filt
# 24583 features across 7044 samples

head(data.filt)
table(data.filt$orig.ident)
DimPlot(data.filt,label = T)

# 加载注释结果
load('inputs/phe.Rdata')
data.filt$celltype <- phe$celltype
DimPlot(data.filt,group.by = 'celltype',label = T)

注释的umap结果如下：

2、过滤1514个基因之后的降为聚类结果

读取基因过滤并进行降维聚类：

# 读取1514个基因
gs <- read.table('remov_gene_list.txt')[,1]
gs

gs_left <- rownames(data.filt)[rownames(data.filt) %in% gs]
gs_left
length(gs_left)
# [1] 232

# 常规方法查看需要过滤的基因
genes_mt =  rownames(data.filt)[grepl("^MT-", rownames(data.filt))]
genes_mt
genes_ribo =  rownames(data.filt)[grepl("^RP[SL]", rownames(data.filt))]
genes_ribo
genes_hsp =  rownames(data.filt)[grepl("^HSP", rownames(data.filt))]
genes_hsp 

data.filt
sce.all.filt <- data.filt[!rownames(data.filt) %in% gs ,]
sce.all.filt
# 24351 features across 7044 samples within 1 assay

###### step2: 基本的降维聚类分群  ###### 
sce <- sce.all.filt %>% 
  NormalizeData( ) %>%  
  FindVariableFeatures( ) %>%  
  ScaleData( ) %>%  
  RunPCA( )  %>% 
  RunHarmony( "orig.ident") %>%
  FindNeighbors(dims = 1:15) %>% 
  FindClusters(resolution = 0.5) %>% 
  RunUMAP(dims = 1:15, reduction = "harmony") 
  
# 对比两次结果
p1 <- DimPlot(data.filt,group.by = 'celltype',label = T) + ggtitle("Before remove")
p2 <- DimPlot(sce ,group.by = 'celltype',label = T) + ggtitle("After remove")
p1+p2

结果如下：基本上没有什么影响。

之所以结果影响非常小，想必是因为降维聚类用到的高变基因里面 这1514个基因的贡献度不大，不知道其他的数据集中结果会怎样！

# 查看数据中 高变基因与这1514个基因的交集
length(VariableFeatures(data.filt))
gs_vsd <- intersect(VariableFeatures(data.filt),gs_left)
gs_vsd
length(gs_vsd)

# 展示其在数据中的变异情况
plot1 <- VariableFeaturePlot(data.filt)
plot2 <- LabelPoints(plot = plot1, points = gs_vsd, repel = TRUE)
plot2

共有43个基因是高变基因，且大部分都在红色的底部：

这43个基因在不同亚群中的表达情况：

# 查看这些基因在不同亚群中的表达情况
DoHeatmap(data.filt, features = gs_vsd, group.by = "celltype")