📚 课程目录第3章真实项目背景

第三章：多批次RNA-seq项目中的注释矛盾

说实话，做转录组项目最怕什么？不是低质量的测序数据，也不是复杂的差异分析。我最怕的是——样本测完，分析做完，结果却对不上。尤其是跨批次的项目，不同时间、不同平台、不同分析人员跑出来的结果，注释信息经常打架。

我个人习惯在项目启动前就问清楚：所有批次的注释文件一致吗？参考基因组版本统一吗？但现实往往很骨感。今天我就跟你聊聊，我在一个真实的癌症转录组项目里，是怎么被注释矛盾坑了一把的。

3.1 问题：基因名对不上号

项目背景很简单：研究某癌种的转录组特征，共100个样本，分三批测序。

第一批：40个样本，由合作方A完成，使用Ensembl GRCh38 v90注释。
第二批：35个样本，自己实验室测的，使用RefSeq GRCh38.p7。
第三批：25个样本，外包给公司B，使用GENCODE v29。

三个批次的原始数据处理流程都不一样。我拿到统一表达矩阵后，发现一个严重问题：同样的基因ID，在不同批次里居然有不同名字。

关键矛盾：Ensembl ID、RefSeq ID、Gene Symbol三者之间的映射存在不一致。比如"BRCA2"在Ensembl里是"ENSG00000139618"，但在RefSeq里对应的是"NM_000059"，而在GENCODE里又变成了"ENST00000380152"这样的转录本ID。

3.2 具体案例：一波三折的基因对比

我记得最清楚的一个基因是TP53。在第一批数据里，表达量用"TP53"这个符号表示，正常。到第二批数据里，同一个基因被标注成了"P53"，第三批更离谱，直接用了"ENSG00000141510"这个Ensembl ID。

如果你不仔细检查，直接把这几个列当成不同基因来算差异表达，结果会怎样？

答案很简单：分析报废。

我当时用了个小脚本，把所有基因的ID映射到统一的Gene Symbol上，才发现了问题。下面是我用来检查注释不一致的Python代码片段：

import pandas as pd

# 读取三个批次的表达矩阵
batch1 = pd.read_csv("batch1_expression.csv", index_col=0)
batch2 = pd.read_csv("batch2_expression.csv", index_col=0)
batch3 = pd.read_csv("batch3_expression.csv", index_col=0)

# 检查基因ID列的命名
print("Batch1 基因名示例:", batch1.index[:5].tolist())
print("Batch2 基因名示例:", batch2.index[:5].tolist())
print("Batch3 基因名示例:", batch3.index[:5].tolist())

# 找出共同基因
common = set(batch1.index) & set(batch2.index) & set(batch3.index)
print(f"共同基因数量: {len(common)}")

# 如果共同基因太少，大概率是注释不一致
if len(common) < 1000:
    print("警告：注释可能不一致！")

跑完以后，共同基因只有482个。正常情况，人类转录组至少应该有15000+个共同蛋白编码基因。你看看，问题严重吧。

3.3 为什么注释会不一致？

原因其实不复杂。每次基因组更新，都会修改一些基因的ID、符号或者坐标。尤其是Ensembl和RefSeq，两个数据库的命名方式完全不同。

我整理了一个表格，把最常见的注释差异列出来：

差异类型	示例	影响范围
Gene Symbol不同	TP53 vs P53 vs ENSG00000141510	跨批次合并失败
染色体坐标偏移	GRCh37 vs GRCh38	reads计数不一致
基因类型定义	蛋白编码 vs lncRNA注释不同	功能注释错误
ID版本号差异	NM_000059 vs NM_000059.4	无法匹配

我这里有个小技巧：在项目开始前，就用统一脚本把所有样本的注释文件转换成同一个格式。我习惯用GENCODE的平级基因注释文件（GTF），因为它同时提供了Gene ID、Gene Symbol和基因类型，映射最稳定。

3.4 解决方案：统一注释源

当时我被迫做了两件事。第一，重新从原始数据开始，统一用Ensembl 105注释文件重新定量所有样本。第二，建立一个从Ensembl ID到最新Gene Symbol的映射表。

流程太长了，我给个示例步骤：

下载同一个版本（Ensembl v105）的GTF文件。
对所有样本的fastq文件，用STAR重跑比对。
用featureCounts统一计数，生成表达矩阵。
用Ensembl自带的ID映射工具（biomaRt）把基因ID转成Gene Symbol。
检查映射后的符号是否还有冲突。

下面是使用R语言的biomaRt包进行ID映射的示例：

library(biomaRt)

# 连接到Ensembl数据库
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 准备要映射的Ensembl ID列表
ensembl_ids <- c("ENSG00000141510", "ENSG00000133703", "ENSG00000156970")

# 映射到Gene Symbol和基因类型
result <- getBM(
    attributes = c("ensembl_gene_id", "hgnc_symbol", "gene_biotype"),
    filters = "ensembl_gene_id",
    values = ensembl_ids,
    mart = ensembl
)
print(result)

你看，这样就能得到统一的命名了。跑完以后，三个批次的共同基因数从482个飙升到18732个。嗯，这才像话。

3.5 批次效应与注释矛盾的叠加效应

你以为统一注释就万事大吉了？别忘了，不同批次的实验条件、湿实验操作、测序深度都不同。这些因素会产生批次效应，叠加在注释差异上，让数据更加混乱。

我曾经做过一次可视化，把三个批次的样本主成分分析（PCA）画在一个图上。结果第一批和第二批的样本完全分离，第三批夹杂在中间。后来发现，除了注释问题，还有测序深度的差异。

核心结论：批次效应和注释矛盾是两种不同的问题，但它们会相互放大。先解决注释矛盾，再处理批次效应，顺序不能搞反。

3.6 一条实用的检查清单

我个人习惯，在合并多批次表达数据前，一定走一遍下面的清单：

第一步：检查所有样本的参考基因组版本是否一致。
第二步：检查所有样本的注释文件（GTF/GFF）是否来自同一数据库和版本。
第三步：统一使用基因ID（优先选择Ensembl Gene ID）作为唯一标识。
第四步：用ID映射工具验证映射后的Gene Symbol是否有一对多的情况。
第五步：做初步的PCA或聚类，看是否有批次相关的结构。

这五步走完，你才能放心地往下做差异分析。

3.7 核心流程示意图

下面这张SVG图概括了注释矛盾的产生与解决流程。我画得比较简洁，但关键步骤都在里面。

嗯，这张图你重点看左侧的“问题来源”和中间的“检查步骤”。只要你走通了从“统一注释源”到“ID映射”再到“检查共同基因数量”这条路径，注释矛盾基本就能避免。

3.8 小结

这次经历给我最大的教训是：不要想当然认为不同批次的注释文件是等价的。我们做数据分析的，最怕的就是数据本身就有问题，而我们还浑然不觉。

注释矛盾听起来是个小问题，但它能直接影响差异表达基因的数量和质量，甚至让你得出完全相反的生物学结论。我建议你在每次分析前，都先跑一遍前面那个检查清单。花半小时排除注释问题，能省掉后面几天的debug时间。

至于更复杂的批次效应移除，比如ComBat、limma、Harmony这些方法，那是下一层的事。先把基础打牢，咱们再谈后面的。