📚 课程目录第4章细胞类型注释

数据格式准备：表达矩阵、元数据与目标基因集

数据格式不对，分析再牛也白搭。

我个人习惯是，拿到单细胞数据后，第一件事不是跑NicheNet，而是检查文件结构。你想想看，如果表达矩阵里基因名是"ENSG000001..."这种ID格式，后面配受体数据库根本匹配不上。我之前见过不少同学上来就跑NicheNet，结果卡在数据读取这一步，一查才发现格式根本不对。

说白了，NicheNet对输入数据有三个硬要求：表达矩阵、元数据（细胞分组/条件）、感兴趣基因集。缺一个都不行。

4.1 表达矩阵：counts还是TPM？

这个问题我常被问到。NicheNet官方推荐用TPM或归一化后的表达值。

为什么？因为配受体分析的核心是比较不同细胞群之间的信号强度。counts数据受测序深度影响太大，同一基因在不同样本里可能差好几倍，这会导致错误的信号推断。当然原始counts数据也不能浪费——计算差异表达基因时它会更有优势。

格式	推荐用途	注意事项
counts	差异表达分析	测序深度影响大，不适合直接用于配受体评分
TPM	配受体分析	基因长度和测序深度已校正，可比性强
log2(TPM+1)	相关性计算	降低高表达基因的权重，适合做相关性分析

避坑：表达矩阵的行名必须是基因名，推荐用HGNC Symbol格式（如"IL6"、"CXCL12"）。如果是Ensembl ID，需要先转换。我曾经因为没注意这个，浪费了一整天的跑图时间。

那么具体怎么准备？我一般这样操作：

# 从Seurat对象中提取表达矩阵
# 假设你已经有了Seurat对象 seurat_obj
# 提取RNA assay的counts数据
counts_mat <- GetAssayData(seurat_obj, assay = "RNA", slot = "counts")
# 转为TPM（这里简化处理，实际项目中建议用更严谨的方法）
library(scater)
tpm_mat <- calculateTPM(counts_mat, 
                         effective_length = gene_lengths)  # 基因长度数据

如果你用的是10X Genomics Cell Ranger的输出，那就更简单了。它直接给了你三个文件：barcodes.tsv、features.tsv、matrix.mtx。用Read10X()读进来后，记得检查一下基因名格式。

# 检查基因名格式
head(rownames(tpm_mat))
# 如果不是HGNC Symbol格式，需要转换
# 我推荐用 biomaRt 包
library(biomaRt)
ensembl <- useEnsembl(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
gene_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
                  filters = "ensembl_gene_id",
                  values = rownames(tpm_mat),
                  mart = ensembl)

4.2 元数据：细胞分组和条件信息

元数据，说白了就是告诉NicheNet哪个细胞属于哪个群。

举个例子，你研究肿瘤微环境，可能有"肿瘤细胞"、"CD8+ T细胞"、"巨噬细胞"这几个群。元数据就是一个表格，第一列是细胞ID，第二列是分组标签。

# 从Seurat对象中提取元数据
metadata <- seurat_obj@meta.data
# 确保包含细胞分组信息
# 比如我们有一个叫 "cell_type" 的列
head(metadata[, c("cell_type", "sample_group")])

# 格式示例
#                 cell_type sample_group
# AAACCCAAGAAACACA    T_cell       tumor
# AAACCCAAGAAACTGA    T_cell       normal
# AAACCCAAGAAACACG Macrophage       tumor

小技巧：如果你的分组信息在另一个文件里，可以用merge()或left_join()合并到Seurat对象的meta.data中。我个人习惯是用dplyr包来处理这种数据合并，代码可读性更好。

元数据还需要包含条件/样本信息。比如你有"疾病组"和"对照组"，或者"处理前"和"处理后"。这些信息在你的数据分析中很关键——NicheNet需要知道哪些细胞群之间有通信关系。

# 添加条件信息
metadata$condition <- ifelse(grepl("tumor", metadata$sample_group), 
                              "disease", "healthy")
# 检查是否所有细胞都有分组
table(metadata$cell_type, metadata$condition)

4.3 感兴趣基因集：从哪里来？

这是最容易被忽视的一步。感兴趣基因集（target gene set），就是你想分析的下游响应基因。

举个例子，你研究T细胞和肿瘤细胞之间的通信。你可能想知道：肿瘤细胞通过哪些配体影响了T细胞？那么T细胞中上调的基因就是你的目标基因集。这些基因可以从FindMarkers()结果中提取。

# 从Seurat中找差异表达基因
# 假设你想找T细胞中肿瘤组 vs 正常组的差异基因
markers <- FindMarkers(seurat_obj,
                       ident.1 = "tumor_T_cell",
                       ident.2 = "normal_T_cell",
                       min.pct = 0.25,
                       logfc.threshold = 0.25)

# 提取显著上调的基因（p_val_adj < 0.05, avg_log2FC > 0.25）
target_genes <- rownames(markers[markers$p_val_adj < 0.05 & 
                                  markers$avg_log2FC > 0.25, ])
# 查看前几个基因
head(target_genes)
# 输出示例：
# [1] "GZMB"  "PRF1"  "IFNG"  "TNF"

注意，目标基因集的质量直接影响分析结果。我建议：

不要只看p值，也要看log2FC，设置合适的阈值
建议目标基因集在20-200个之间，太少没法做富集分析，太多会引入噪声
如果你有先验知识（比如某个通路的关键基因），可以手动补充进去

    实战经验：我在一个项目里同时用了三个不同阈值（log2FC>0.25、>0.5、>1.0）的目标基因集做NicheNet分析，结果发现中等阈值（>0.5）的预测结果与后续实验验证最吻合。这里没有绝对标准，建议你多试试。

4.4 数据准备的整体流程

下面这张图可以帮你理清脉络。简单说就是三路数据汇合，统一格式后输入NicheNet。

4.5 常见问题与避坑指南

讲到这里，我提几个实际项目中踩过的坑：

细胞ID不匹配：表达矩阵的列名和元数据的行名必须完全一致。一个常见错误是元数据里的细胞ID多了引号或空格。
重复基因名：有些基因有多个别名，比如"TNF"和"TNFA"其实是一个基因。我建议用unique()检查一下，重复的基因取平均表达值。
零表达基因：如果某个基因在所有细胞中都不表达，留着也没意义。建议过滤掉全零行。

# 检查并修复常见问题
# 1. 检查重复基因名
sum(duplicated(rownames(expr_mat)))
# 如果有重复，用 aggregate 或 tapply 合并

# 2. 过滤全零基因
expr_mat <- expr_mat[rowSums(expr_mat) > 0, ]

# 3. 确保细胞ID匹配
stopifnot(all(colnames(expr_mat) %in% rownames(metadata)))
# 如果匹配不上，调整列的顺序
expr_mat <- expr_mat[, rownames(metadata)]

注意：千万不要在NicheNet分析前对表达矩阵做log归一化，它会默认你的输入是RNA counts或TPM。做过log转换的数据会让配受体分数计算完全偏离预期。

嗯，到这里数据准备就差不多了。说白了，这个环节看着简单，但做不好后面全是白费功夫。我之前有个学生就是元数据里细胞ID多了一个下划线，愣是找了两天bug才排查出来。你按我说的这几步走，基本不会出问题。