单细胞转录组数据挖掘：差异基因分析与功能富集实战

在单细胞转录组测序（scRNA-seq）技术日益成熟的今天，如何从海量数据中挖掘有价值的信息，找到细胞类型特异的差异基因，并理解这些基因的功能，成为了生物信息学分析的关键环节。本文将深入探讨差异基因分析和富集分析在单细胞转录组研究中的应用，并提供实战案例和避坑经验。

差异基因分析：寻找细胞间的秘密

什么是差异基因？

差异基因（Differentially Expressed Genes, DEGs）是指在不同细胞群体或实验条件下表达水平显著不同的基因。在单细胞转录组分析中，我们的目标是识别在不同细胞类型或处理组之间表达量存在显著差异的基因。这些基因往往与特定细胞的功能、发育过程或疾病状态密切相关。

差异基因分析的常用方法

常用的差异基因分析方法包括：

• t 检验和方差分析 (ANOVA)：经典的统计学方法，适用于比较两组或多组样本的均值差异。但需要注意，单细胞数据通常不满足正态分布的假设，因此需要进行数据转换或使用非参数检验。
• Wilcoxon 秩和检验 (Mann-Whitney U test)：一种非参数检验，适用于比较两组非正态分布的数据。在单细胞分析中，常用于比较两组细胞类型的基因表达差异。
• DESeq2 和 edgeR：最初为 bulk RNA-seq 数据设计的差异表达分析工具，经过改进后也可用于单细胞数据。它们基于负二项分布模型，能够有效地处理 count 数据，并考虑样本间的文库大小差异。
• limma-voom：另一种为 RNA-seq 数据设计的差异表达分析工具，通过 voom 函数将 count 数据转换为 log-counts per million (log-CPM) 值，并使用线性模型进行分析。也适用于单细胞数据。
• scran 和 MAST：专门为单细胞数据设计的差异表达分析工具。scran 提供了 normalization 的方法，MAST (Model-based Analysis of Single-cell Transcriptomics) 考虑了单细胞数据中常见的 dropout 事件，并使用 hurdle 模型进行分析。

实战代码：使用 Seurat 进行差异基因分析

Seurat 是一个常用的单细胞分析 R 包。以下代码展示了如何使用 Seurat 查找两个细胞群之间的差异基因：

# 安装 Seurat (如果尚未安装)
# install.packages('Seurat')

library(Seurat)

# 假设你已经创建了一个 Seurat 对象，名为 'seurat_object'
# 例如，可以使用 Read10X 或 ReadH5AD 函数读取数据

# seurat_object <- Read10X(data.dir = "path/to/your/data")
# seurat_object <- CreateSeuratObject(counts = seurat_object, project = "scRNAseq")

# 假设你已经进行了数据质控、normalization 和聚类分析

# 查找细胞类型 1 和细胞类型 2 之间的差异基因
# 使用 Wilcoxon 秩和检验
degs <- FindMarkers(seurat_object, ident.1 = "cell_type_1", ident.2 = "cell_type_2", test.use = "wilcox")

# 查看结果
head(degs)

# 可以设置 logfc.threshold 参数来过滤显著差异的基因
degs_significant <- FindMarkers(seurat_object, ident.1 = "cell_type_1", ident.2 = "cell_type_2", test.use = "wilcox", logfc.threshold = 0.25)
head(degs_significant)

#可以将结果保存到csv文件
write.csv(degs,file="diff_genes.csv")

避坑指南

• 数据质控至关重要：在进行差异基因分析之前，必须对原始数据进行严格的质控，包括去除低质量细胞、过滤低表达基因等。这会直接影响后续分析的准确性。
• Normalization 方法的选择：不同的 normalization 方法适用于不同的数据集。常见的 normalization 方法包括 CPM, TPM, RPKM, 以及 Seurat 中常用的 LogNormalize。需要根据数据的特点选择合适的 normalization 方法。
• 多重检验校正：由于差异基因分析涉及大量的假设检验，需要进行多重检验校正，以控制假阳性率。常用的校正方法包括 Bonferroni 校正和 Benjamini-Hochberg (FDR) 校正。
• 细胞类型注释的准确性：差异基因分析结果的解读依赖于细胞类型注释的准确性。建议使用多种方法进行细胞类型注释，并结合已知的生物学知识进行验证。

富集分析：揭示基因背后的功能

什么是富集分析？

富集分析（Enrichment Analysis）是一种用于确定一组基因或蛋白质是否在特定的生物学通路或功能类别中过度表示的方法。通过富集分析，我们可以了解差异基因所参与的生物学过程，从而更好地理解细胞的功能和调控机制。

常用的富集分析方法

常见的富集分析方法包括：

• GO 富集分析 (Gene Ontology Enrichment Analysis)：GO (Gene Ontology) 是一种对基因和蛋白质功能进行标准化的分类体系。GO 富集分析用于确定差异基因是否在特定的 GO term 中过度表示。
• KEGG 富集分析 (Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis)：KEGG 是一种对生物学通路进行分类的数据库。KEGG 富集分析用于确定差异基因是否在特定的 KEGG pathway 中过度表示。
• GSEA (Gene Set Enrichment Analysis)：GSEA 是一种基于基因集的富集分析方法。它不需要预先筛选差异基因，而是直接使用所有基因的表达数据进行分析。

实战代码：使用 clusterProfiler 进行富集分析

clusterProfiler 是一个常用的 R 包，用于进行 GO 和 KEGG 富集分析。以下代码展示了如何使用 clusterProfiler 对差异基因进行富集分析：

# 安装 clusterProfiler (如果尚未安装)
# if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
# BiocManager::install("clusterProfiler")

library(clusterProfiler)
library(org.Mm.eg.db) # Mouse 数据库。Human 数据库使用 org.Hs.eg.db

# 假设你已经获得了一组差异基因的基因 ID，存储在 'gene_list' 向量中
# gene_list <- rownames(degs_significant)

# 将基因名转换为 Entrez ID
# 需要根据物种选择合适的数据库，例如 org.Hs.eg.db (Human) 或 org.Mm.eg.db (Mouse)

gene_list <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)

# 进行 GO 富集分析
go_enrichment <- enrichGO(gene          = gene_list$ENTREZID,
                        OrgDb         = org.Mm.eg.db,
                        keyType       = 'ENTREZID',
                        ont           = "BP", # Biological Process
                        pAdjustMethod = "BH",
                        pvalueCutoff  = 0.05,
                        qvalueCutoff  = 0.05,
                        readable      = TRUE)

# 查看结果
head(go_enrichment)

# 进行 KEGG 富集分析
kegg_enrichment <- enrichKEGG(gene          = gene_list$ENTREZID,
                          organism      = 'mmu', # Mouse
                          keyType       = 'kegg',
                          pvalueCutoff  = 0.05,
                          pAdjustMethod = "BH",
                          qvalueCutoff  = 0.05)

# 查看结果
head(kegg_enrichment)

# 可以使用 dotplot 或 barplot 函数可视化富集分析结果
dotplot(go_enrichment, showCategory = 10)
barplot(kegg_enrichment, showCategory = 10)

避坑指南

• 选择合适的数据库：GO 和 KEGG 数据库会不断更新，需要选择最新版本的数据库进行分析。另外，需要根据物种选择合适的数据库，例如 Human 对应 org.Hs.eg.db，Mouse 对应 org.Mm.eg.db。
• ID 转换的重要性：富集分析需要使用标准的基因 ID，例如 Entrez ID。如果你的基因列表使用的是其他类型的 ID，需要先进行 ID 转换。
• p 值和 FDR 的理解：富集分析会计算每个 term 的 p 值，用于衡量该 term 是否显著富集。由于涉及大量的假设检验，需要进行多重检验校正，得到 FDR 值 (False Discovery Rate)。通常使用 p 值 < 0.05 或 FDR < 0.05 作为显著性阈值。
• 结果的生物学意义：富集分析的结果需要结合已知的生物学知识进行解读。需要仔细分析富集到的 term，判断它们是否与你的研究问题相关。

总结

差异基因分析和富集分析是单细胞转录组数据分析的重要步骤。通过这两个步骤，我们可以识别细胞类型特异的基因，并理解这些基因的功能。在实际应用中，需要注意数据质控、Normalization 方法的选择、多重检验校正以及生物学意义的解读。掌握这些技巧，将能更有效地利用单细胞转录组数据，探索生命科学的奥秘。