公共数据常见坑与避雷

跟完这篇，您能在拿到任何一份公共数据后，按顺序排掉四个最容易翻车的坑——单位喂错、样本贴错、探针没注释、批次没识别——确认矩阵干净、分组对得上、能直接开跑，而不是分析做到一半才发现数据从头就错了。

先把这四个坑记在心上

公共数据（GEO、TCGA、各类补充文件）省了您测序的钱，但也把一堆隐患一起打包给了您。它们的共同点是：不报错、看着正常，结果却是错的。下面这张表把四个坑、它们的表现和后果摆在一起，本篇后面每一节各拆一个。

坑	典型表现	不查的后果
TPM 当成 counts	矩阵是小数、每列加起来≈1e6	差异分析报错，或悄悄算错
列名 ≠ 样本顺序	矩阵和分组表分开下载	分组贴反，差异基因全是假的
探针没注释	行名是 `1007_s_at` 这类探针号	一整张结果全是探针、没法解读
批次没识别	数据来自多平台、多批次、多中心	把技术差异当成了生物学差异

建议把这四步当成拿到数据后的排雷自查，按顺序过一遍再动分析。

坑一：TPM 不是 counts

公共数据里投稿者上传的，很多是已经归一化的 TPM / FPKM，而不是原始计数。基于 count 的差异分析（DESeq2、edgeR）必须要原始整数 counts，它们内部自带归一化，喂进 TPM 要么直接报错，要么不报错却把方差估计全带偏——后者更危险。

拿到一份矩阵，先花十秒判断它是什么单位，看两个特征就够：

# mat：读进来的表达矩阵
all(mat == round(mat))   # 全是整数 → 大概率是原始 counts
head(colSums(mat))       # 每列都≈1e6 → 是 TPM；各列不等且带小数 → FPKM / CPM

全是非负整数 → 基本是 raw counts，可以直接做差异分析。
带小数、且每列和都约等于 1e6 → 是 TPM，只能拿来看和展示，不能当差异分析的输入。
带小数、每列和各不相同 → 多半是 FPKM 或 CPM，同样别喂给 DESeq2。

拿不准或只拿到归一化值时，最稳的做法是回上游要一份原始 counts（featureCounts、HTSeq、Salmon/kallisto 的 counts 输出都行）。四种单位分别归一化了什么、各自用在哪，见表达量单位怎么选。

坑二：矩阵列名和样本顺序对不上

这是最隐蔽、后果最重的坑。公共数据里，表达矩阵和样本信息常常是两个文件分开下载的：矩阵在补充文件里，分组 / 表型在 series_matrix 或另一张临床表里。两者的样本排列顺序不保证一致。

新手最容易犯的错，是默认「矩阵的第 i 列 = 分组表的第 i 行」，直接把分组标签贴上去：

# 危险写法：假设两边顺序天然一致，直接拼
coldata <- data.frame(condition = meta$condition)   # 顺序错了，分组就贴反了

顺序一旦错位，tumor 被当成 normal，跑出来的差异基因看着漂亮、其实全是假的，而且全程不报任何错。

正确做法是永远按样本名对齐，而不是按位置：先确认两边样本集合一致，再按分组表的顺序重排矩阵的列。

# 先确认两边是同一批样本，一个不多一个不少
setdiff(colnames(mat), meta$sample)   # 应为 character(0)
setdiff(meta$sample, colnames(mat))   # 应为 character(0)

mat <- mat[, meta$sample]              # 按分组表顺序重排矩阵的列
all(colnames(mat) == meta$sample)      # 必须为 TRUE，才能往下

两个 setdiff 都返回 character(0)、最后一行返回 TRUE，才算真正对齐。名字对不上时，多半是大小写、分隔符（- 和 . 和 _）、或末尾多了空格——读文件时加 check.names = FALSE，能避免 R 偷偷把 Sample-1 改成 Sample.1。对齐的完整做法见表达矩阵与分组表整理。

坑三：芯片探针没注释

如果数据来自表达芯片（microarray），矩阵的行名往往不是基因，而是探针号，像 1007_s_at、ILMN_1651209。您得先借助平台（GPL）的注释表，把探针映射成基因 symbol，否则后面画出来的火山图、热图，坐标轴上全是没人看得懂的探针号。

第一步，取到对应平台的注释表。用 GEOquery 按 GPL 编号拉：

library(GEOquery)
gpl  <- getGEO("GPL570")          # 换成您数据实际用的平台编号
anno <- Table(gpl)                # 探针号与基因的对照表
colnames(anno)                    # 先看清 symbol 在哪一列，各平台列名不一样

注意：注释表里放 symbol 的列，各平台命名不统一（有的叫 Gene Symbol，有的叫 GENE_SYMBOL 或 Symbol），所以先 colnames() 看清楚再取，别照抄别人的列名。除了 GPL 表，也可以用对应平台的 Bioconductor 注释包。

第二步，把探针换成基因。同一个基因常对应多个探针，行名又不能重复，得去重。一个被广泛采用的惯例是保留平均表达最高的那个探针（是常用做法、按需也可换成取中位数或用专门的合并函数，不是铁律）：

sym <- anno[["Gene Symbol"]][match(rownames(mat), anno$ID)]   # 探针 → symbol
keep <- !is.na(sym) & sym != ""                               # 丢掉没注释上的探针
mat  <- mat[keep, ]; sym <- sym[keep]

ord  <- order(-rowMeans(mat))          # 按平均表达从高到低排
mat  <- mat[ord, ]; sym <- sym[ord]
mat  <- mat[!duplicated(sym[ord]), ]   # 每个基因只留表达最高的那一行
rownames(mat) <- sym[ord][!duplicated(sym[ord])]

RNA-seq 也有一个近亲坑：行名是带版本号的 Ensembl ID（ENSG00000141510.12），做 ID 转换或跨数据集合并前，先把版本后缀去掉，否则同一个基因会因为版本号不同而对不上：

rownames(counts) <- sub("\\..*$", "", rownames(counts))   # 去掉 .12 这类版本后缀

坑四：批次没识别

公共数据常常把不同平台、不同实验室、不同时间测的样本混在一起。批次（batch）带来的技术差异，可能比您关心的生物学差异还大。如果不识别它，您以为找到的「肿瘤 vs 正常」差异，可能只是「A 批次 vs B 批次」的机器差异。

识别批次最直接的办法是画 PCA 看样本怎么抱团：按批次上色，如果样本是按批次分开、而不是按分组分开，批次效应就压过了生物学信号。

library(DESeq2)
vsd <- vst(dds)                                      # 方差稳定变换后再看结构
plotPCA(vsd, intgroup = c("batch", "condition"))     # 分别按批次、按分组上色对比

看到批次效应后，怎么处理分两种情况：

情况	判断	怎么办
批次与分组没有混杂	每个批次里两组都有	把批次当协变量放进模型
批次与分组完全混杂	一个批次全是肿瘤、另一个全是正常	无解，任何方法都救不回来

第一种情况，把批次写进差异分析的设计公式，让模型在估计分组效应时顺带扣掉批次——这是标准做法：

design(dds) <- ~ batch + condition   # 把 batch 当协变量，和分组一起建模
dds <- DESeq(dds)

如果只是想画一张干净的 PCA 或热图，可以用 limma 的 removeBatchEffect 得到去批次后的矩阵——但它只用于可视化和聚类，去批次后的值不能拿去做差异分析（差异分析要走上面把 batch 放进设计公式的路子）：

library(limma)
logexpr_adj <- removeBatchEffect(logexpr, batch = batch)   # 仅供画图 / 聚类，别当差异分析输入

诚实说一句边界：如果批次和分组完全混杂（比如所有肿瘤样本都在一个批次、所有正常样本都在另一个批次），那么批次效应和生物学效应在数学上分不开，没有任何统计方法能挽救——这种数据集只能放弃或换一批。PCA、火山、热图这些图怎么读，见结果怎么读。

上手前的四步自查清单

把上面四节收成一张可以照着打勾的单子，拿到任何公共数据，动分析前先过一遍：

单位：all(mat == round(mat)) 是不是 TRUE、列和是不是≈1e6，确认是原始 counts 还是归一化值。
对齐：两个 setdiff 都返回 character(0)、all(colnames(mat) == meta$sample) 为 TRUE，按名字而非位置对齐。
注释：行名是基因还是探针；是探针就映射成 symbol，是带版本号的 Ensembl 就去掉后缀。
批次：画一张 PCA，按批次上色看有没有抱团；有就查是否与分组混杂，再决定放进模型还是放弃。

这四步都打勾，数据才算「干净」，可以放心接差异分析。

在平台上跑一遍

在百沐一下，这几步排雷不用您写代码：把矩阵和样本信息传上去，组学智能体会自动判断单位是不是原始 counts、按样本名对齐两张表、认出芯片探针并映射成基因、再用 PCA 帮您看有没有批次抱团，把可疑的地方标出来请您确认。上传数据、说清目的即可：组学智能体使用说明，或直接开始 app.baimuyixia.com。