跳到主要内容

03 数据探索：拿到新数据先看什么

数据到手了，最大的诱惑是直接跑差异分析。但这是新手最容易翻车的地方 — 不探索就分析，等于闭着眼开车。

本章的目标：花 1-2 小时做一轮探索性分析，回答三个问题：

  1. 数据质量行不行？ — 有没有明显的离群样本或批次效应
  2. 临床变量分布怎么样？ — 样本够不够做我想做的比较
  3. 初步信号在哪？ — 主成分和哪些临床变量相关

这三个问题答完再决定分析路线，能避免 80% 的"做到一半发现数据有问题"的悲剧。

加载上一章保存的数据​

library(DESeq2)
library(ggplot2)
library(dplyr)

# 上一章保存的三个文件
se <- readRDS("data/tcga_lihc_rnaseq.rds")
clin <- read.csv("data/tcga_lihc_clinical.csv")

第一步：看表达矩阵的整体分布​

# 只取 tumor 样本（barcode 14-15 位 = "01"）
tumor_idx <- substr(colnames(se), 14, 15) == "01"
se_tumor <- se[, tumor_idx]

# 构建 DESeq 对象（这里只做 vst，不做差异）
dds <- DESeqDataSet(se_tumor, design = ~ 1)
keep <- rowSums(counts(dds) >= 10) >= 10
dds <- dds[keep, ]

# VST 变换
vsd <- vst(dds, blind = TRUE)

库大小分布​

lib_sizes <- colSums(counts(dds)) / 1e6

ggplot(data.frame(size = lib_sizes), aes(x = size)) +
  geom_histogram(bins = 30, fill = "steelblue") +
  labs(x = "Library size (millions)", y = "Samples",
       title = "TCGA-LIHC: Library size distribution") +
  geom_vline(xintercept = median(lib_sizes), linetype = "dashed", color = "red")

看什么：大部分样本应该在同一个量级（比如都在 20-80M）。如果有样本 <5M 或 >200M，标记为异常。

检出基因数​

genes_per_sample <- colSums(counts(dds) > 0)
summary(genes_per_sample)

第二步：PCA — 全局质量检查​

PCA 是探索阶段的核心工具。不是为了找差异基因，而是为了回答"样本之间的主要差异来源是什么"。

pca_data <- plotPCA(vsd, intgroup = "1", returnData = TRUE)

# 把临床信息合并进来
pca_data$patient_id <- substr(rownames(pca_data), 1, 12)
pca_merged <- left_join(pca_data, clin, by = "patient_id")

# 按分期上色
ggplot(pca_merged, aes(x = PC1, y = PC2, color = stage)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "PCA: colored by stage") +
  theme_minimal()

PCA 上要检查的事​

看什么正常异常
有没有离群点所有样本在一个主群里某个点远离群体 → 检查这个样本的 QC
PC1 和什么相关和你关心的生物学变量相关和你不关心的变量（批次、性别）相关 → 需要校正
分组是否分开你要比较的两组在 PC1/2 上有分离完全混在一起 → 效应可能很弱
# PC1 vs 每个临床变量的相关
cor.test(pca_merged$PC1, pca_merged$age)
wilcox.test(pca_merged$PC1 ~ pca_merged$gender)
kruskal.test(pca_merged$PC1 ~ pca_merged$stage)

如果 PC1 和 batch / plate 显著相关，说明批次效应是数据里最大的变异来源，必须在差异分析时控制。

第三步：临床变量分布​

不看临床分布就做分析，会遇到"治疗组全是男性、对照组全是女性"这种混杂。

# 分期分布
table(clin$stage, useNA = "ifany")

# 性别 × 分期交叉表
table(clin$gender, clin$stage)

# 生存时间分布
ggplot(clin, aes(x = OS.time / 30, fill = factor(OS))) +
  geom_histogram(bins = 30, position = "stack") +
  labs(x = "Follow-up (months)", y = "Patients",
       fill = "Event", title = "OS time distribution")

要重点检查的​

检查项为什么
每组样本数少于 20 的组做差异分析没统计力
事件数（死亡数）KM 生存分析至少要 20 个事件
混杂因素分期、性别、年龄在你的分组里是否平衡
缺失比例关键变量缺失 >30% 要考虑是否还能用

第四步：表达特征的初步探索​

Top 变异基因​

rv <- rowVars(assay(vsd))
top_var <- order(rv, decreasing = TRUE)[1:500]

# 看 top 500 变异基因里有没有已知的 marker
library(org.Hs.eg.db)
top_symbols <- mapIds(org.Hs.eg.db,
                      keys = rownames(vsd)[top_var],
                      keytype = "ENSEMBL", column = "SYMBOL")
head(sort(top_symbols), 20)

如果 top 变异基因全是线粒体 / 核糖体 / 免疫球蛋白，说明这批数据的主要变异来源可能不是肿瘤本身，而是样本组成差异。

已知 marker 检查​

# 肝癌经典 marker
markers <- c("GPC3", "AFP", "ALB", "EPCAM", "KRT19")
marker_ens <- mapIds(org.Hs.eg.db, keys = markers,
                     keytype = "SYMBOL", column = "ENSEMBL")

# 在 tumor vs normal 里看这些 marker
boxplot(assay(vsd)[marker_ens["GPC3"], tumor_idx],
        assay(vsd)[marker_ens["GPC3"], !tumor_idx],
        names = c("Tumor", "Normal"), main = "GPC3 expression")

GPC3 在肿瘤里应该显著高于 normal。如果不是，检查数据或样本标签。

探索完成后做一个判断​

花了 1-2 小时，现在应该能回答：

问题回答示例
数据质量"有 2 个离群样本需要删除，其余正常"
批次"PC1 和 plate 无关，不需要额外校正"
样本量"Stage III+IV 只有 90 人、60 个事件，做生存够用"
初步信号"PC1 和 stage 弱相关（p=0.03），说明分期信号存在但不强"
下一步"可以按计划做差异 + 生存 + 免疫浸润"

如果探索结果和预期严重不符（比如 PCA 完全没分期信号），要回头考虑是选题需要调整还是数据有问题。

常见坑​

坑 1：不做 VST 直接画 PCA​

Raw counts 上画 PCA 没有意义 — 会被高表达基因主导。VST 或 rlog 是必须的前处理步骤。

坑 2：把探索当分析，p 值直接报​

探索阶段的 p 值都是"看一眼方向用的"，不能直接报在论文里。要报告的 p 值必须来自正式的差异分析（有设计、有校正）。

坑 3：发现离群样本直接删​

先搞清楚为什么离群：库大小太低？样本标签错了？组织类型不对？有合理理由才能删，不能"看着不顺眼就删"。方法段要写清删除原因。

坑 4：不看事件数就做生存分析​

371 个 tumor 听上去样本量够，但如果只有 30 个人死亡（事件数 = 30），KM 曲线的置信区间会很宽。TCGA-LIHC 大约有 130 个死亡事件，做二分组生存分析够用，做更精细的分层会吃力。

坑 5：临床变量有混杂不处理​

发现 Stage III 患者大多是男性、Stage I 大多是女性（假设）。这时候"Stage vs 预后"的分析会被性别混杂。处理方法：多变量 Cox 回归加 gender 为协变量，或在单变量分析时注明 limitation。

下一步​

接着深入：

横向延伸：

参考资源​

AI 陪学

让 AI 陪我学这一篇

AI 会读这篇文章后给你 3-5 步学习计划， 逐步陪你学完，最后出 1-3 道题验证你掌握得怎么样。 登录后 AI 才能记住你的进度。

静态文件

离线资料下载

手册 HTML / PDF 已在后台预生成，点击后直接下载网站静态资源。

本篇 HTML本篇 PDF说明