# 第 06 章 数据获取

> 本章属于“Track B：数据基础”。环境搭好（第 05 章）之后，我们正式开始下载本课题的两套核心数据：**TCGA-LIHC**（发现队列）与 **GSE14520**（验证队列）。

## 本章目标

学完本章，你将能够：

1. 从 UCSC Xena 网页下载 TCGA-LIHC 的表达矩阵与临床表（并了解 GDC 的下载思路）；
2. 用 R 的 GEOquery 包下载 GSE14520，提取表达矩阵与临床信息；
3. 理解 GSM→GSE 的层级关系，以及探针 ID 如何映射到基因 Symbol；
4. 按规范组织 `data/` 目录（原始数据只读、脚本与结果分开）；
5. 对下载的数据做最简单的质控初检（维度、缺失值、分组列）。

> 重要提示：本章所有下载都需要**联网**，且 TCGA/GEO 文件较大，下载可能较慢；教学演示时也可以只取部分样本（见 6.2 节提示）。

## 6.1 下载 TCGA-LIHC：推荐用 UCSC Xena

开始前请先做两件事：一是在 `data/raw/` 下建好 `TCGA-LIHC/` 与 `GSE14520/` 两个文件夹（见 6.3 节）；二是确认网络稳定——TCGA 表达矩阵文件较大，中途断网可能导致下载不完整，建议在网速较好的时段下载。

TCGA 数据可以从官方 GDC 下载，但对新手来说，**UCSC Xena**（https://xena.ucsc.edu ）更友好：它把 TCGA 整理成“基因 × 样本”的表达矩阵和“样本 × 字段”的临床表，点几下就能下载。

**UCSC Xena 网页下载步骤**（以当前界面为例；网站界面会更新，以官网为准）：

1. 打开 https://xena.ucsc.edu ，点击 **Launch Xena** 进入数据浏览界面（Xena Browser）；
2. 在左侧数据列表中找到 **TCGA Liver Cancer (LIHC)**，点击展开；
3. 勾选表达数据：选择 **HTSeq - FPKM**（FPKM 是一种按基因长度和测序深度标准化的表达量，适合基因间比较），这是本教程后续分析使用的表达矩阵；如果需要做差异表达用的原始计数，可另下载 **HTSeq - Counts**；
4. 勾选临床数据：选择 **phenotype**（临床表），里面包含生存状态、生存时间、TNM 分期等字段；
5. 使用页面的 **Download** 功能分别下载表达矩阵与临床表。

**下载到的是什么文件？**

- **表达矩阵**：一个压缩的 TSV 文本文件（.gz），第一列是基因 ID（如 ENSG000001…… 或基因 Symbol），之后每一列是一个样本。样本名形如 `TCGA-BC-A10Q-01A`，其中 `-01` 一般代表肿瘤组织、`-11` 代表癌旁正常组织（TCGA 样本类型编码规则，以官方说明为准）；
- **临床表**：样本 × 字段的表格，每个样本一行，包含生存和临床信息。

**FPKM 还是 Counts？** 本教程选用 HTSeq-FPKM：它按基因长度与测序深度做了标准化，适合比较基因之间的表达水平，也是 limma 差异表达分析（第 07 章）的常见输入（通常先做 log2 转换）。HTSeq-Counts 是原始计数，更适合 DESeq2 / edgeR 等专门工具，本教程暂不涉及。

**GDC 简述（了解即可）**：https://portal.gdc.cancer.gov → Repository → 左侧筛选（Cancer Type 选 Liver Hepatocellular Carcinoma；Data Category 选 Transcriptome Profiling；Experimental Strategy 选 RNA-Seq）→ 把符合条件的文件加入购物车 → 下载表达矩阵与临床文件。GDC 更权威，但筛选项多、文件格式多样，新手阶段先用 Xena。

**下载后检查**：表达矩阵文件通常有几十 MB（解压后更大），临床表只有几百 KB。下载完成后建议先确认文件能正常打开：用记事本或 Excel 打开前几行看看格式（注意：大文件用 Excel 打开会很慢，可用 R 的 `read.table` 或 `data.table::fread` 读取）。文件命名建议加上日期，例如 `LIHC_HTSeqFPKM_2025-01-15.txt.gz`。如果只是想练习完整流程，也可以先用较小的临床表文件熟悉格式。

## 6.2 用 R 的 GEOquery 下载 GSE14520

GSE14520 存放在 GEO 中。GEO 提供了官方 R 包 **GEOquery**，一条命令即可下载整个数据集。

**代码 6-1：下载并读取 GSE14520**

```r
# 第 05 章已安装：BiocManager::install("GEOquery")
library(GEOquery)

# 下载 GSE14520 的 Series Matrix 文件（需联网，视网速可能需要数分钟）
gse <- getGEO("GSE14520", GSEMatrix = TRUE)

# gse 是一个列表：一个 GPL 平台对应一个元素
length(gse)   # 一般等于该 GSE 使用的平台数量
gse[[1]]      # 第一个平台的 ExpressionSet 对象

# 提取表达矩阵：行是探针 ID，列是样本（GSM 编号）
expr <- exprs(gse[[1]])
dim(expr)     # 查看维度，例如 2 万余个探针 × 数百个样本

# 提取样本信息（临床 / 分组 / 随访）
pdata <- pData(gse[[1]])
colnames(pdata)  # 先看有哪些列，再决定怎么用
head(pdata)      # 查看前几行
```

**ExpressionSet 是什么？** GEOquery 把每个平台的完整数据封装成一个 ExpressionSet 对象，它像一个“文件夹”，里面装了三样东西：表达矩阵（用 `exprs()` 取出）、样本信息（用 `pData()` 取出）、探针注释（用 `fData()` 取出）。理解这个结构后，GEOquery 的日常用法就只剩这三个函数了。

**解释几个关键点：**

- `getGEO("GSE14520", GSEMatrix = TRUE)` 下载的是 GEO 官方整理的 Series Matrix 文件；返回的 `gse` 是**列表**，每个元素对应一个平台（GPL）的 ExpressionSet 对象。同一个 GSE 如果用了不同芯片平台，就会有多个元素，所以我们用 `gse[[1]]` 取第一个；
- `exprs()` 取出表达矩阵（ExpressionSet 的“表达”槽），`pData()` 取出样本信息（phenotype data 槽）——这是 GEOquery 最常用的两个函数；
- **GSM→GSE 层级回顾**（第 03 章）：GSE14520 是一个“系列”，包含成百上千个 GSM 样本；`expr` 的列名就是 GSM 编号，`pdata` 的每一行对应一个 GSM；
- 样本分组信息通常在 `pdata` 的 `characteristics_ch1` 等列中，例如 `tissue: tumor` / `tissue: normal`；GSE14520 还带有生存随访字段（这是它能做生存验证的前提）。**不同数据集的列名不同，务必先用 `colnames(pdata)` 查看实际列名**，再写代码提取。

> 小提示：getGEO 会把下载的 Series Matrix 文件缓存在当前工作目录，重复运行时直接读取缓存，速度更快。如果下载中断，删除缓存文件后重新运行即可。

**教学演示提示**：GSE14520 样本数较多，教学时可以先只取部分样本跑通流程：

```r
# 示例：只保留前 20 个样本（教学演示用；正式分析请用全部样本）
expr_sub <- expr[, 1:20]
```

**关于临床列**：GSE14520 的 `pData` 中，样本分组列（tumor / normal）和生存列（生存时间、生存状态）的列名在不同版本中可能不同，常见命名有 `characteristics_ch1`、survival time、survival status 等。第 11 章做生存分析前，我们会先把这些列整理成统一格式：时间列必须是数值（如 12.5 个月），状态列必须是 0/1（0 = 删失 censored，1 = 事件发生）。如果状态列是文字（如 alive / dead），需要先转换为 0/1。

**探针 ID → 基因 Symbol 的映射（平台注释）**

GSE14520 使用的芯片平台是 **GPL3921**（Affymetrix 人类全基因组芯片，具体型号以 GEO 页面为准），表达矩阵的行名是探针 ID（形如 `1552777_at`），不是基因名。要得到基因 Symbol，需要平台注释。三种常用做法：

方法 A（最简单）：Series Matrix 文件本身常带注释列，直接查看：

```r
fData(gse[[1]])                # 探针注释表
colnames(fData(gse[[1]]))      # 看看有哪些注释列（如 Gene Symbol）
```

方法 B（从 GEO 下载 GPL 注释表）：

```r
gpl <- getGEO("GPL3921")       # 下载平台注释
colnames(Table(gpl))           # 查看注释表的列，通常有 ID、Gene Symbol 等
head(Table(gpl))               # 查看前几行示例
```

方法 C（Bioconductor 注释包）：GPL3921 这类 Affymetrix 平台一般有对应的注释包（如 hthgu133a.db，具体以 Bioconductor 当前文档为准），可用 `mapIds()` 做探针→Symbol 的映射。第 07 章预处理时会给出完整映射代码，这里先理解原理：**探针 ID 必须靠平台注释变成基因 Symbol，分析才有生物学意义**。

无论用哪种方法，最终目的都是得到一张“探针 ID → 基因 Symbol”的对照表，然后把它合并到表达矩阵上，把同一基因对应的多条探针合并（通常取平均值或最大值，第 07 章会给出处理重复基因的完整代码）。这一步是芯片数据分析的关键转折点：**从此以后，你的分析对象从“探针”变成了“基因”**，后面所有分析（差异表达、富集、生存）都在基因层面进行。

## 6.3 数据文件怎么组织

下载到的原始文件不要乱放，按第 05 章的约定组织（这是本教程所有章节的统一约定）：

```
hcc_pyroptosis/
├── data/
│   ├── raw/                    # 原始下载文件（只读，永不手动修改）
│   │   ├── TCGA-LIHC/          # Xena 下载的表达矩阵与临床表
│   │   └── GSE14520/           # GEO 下载的系列矩阵文件
│   └── processed/              # 处理后的数据（第 07 章输出）
├── scripts/                    # 01_download.R、02_dea.R ……
├── results/                    # 分析结果表
└── figures/                    # 图片
```

把 Xena 下载的文件放进 `data/raw/TCGA-LIHC/`；把 GEOquery 下载的对象保存到本地，避免以后反复联网下载。`processed/` 目录则留给第 07 章的预处理产物（如探针映射为基因 Symbol、去重、log2 转换后的表达矩阵）——**原始文件永远只读，一切处理都在脚本里可重现**：

```r
# 保存 GSE14520 原始对象到本地（约几十 MB）
save(gse, file = "data/raw/GSE14520/gse14520_raw.RData")
```

**为什么原始数据必须只读？** 分析的可信度来自“可追溯”：任何处理都应当通过脚本完成，而不是手动修改原始文件。一旦你在 Excel 里手动改过原始矩阵，别人（包括半年后的你自己）就无法知道你改了什么，结果也就无法复现。记住一句话：**原始数据只读，脚本重现一切**。

## 6.4 数据质控初检

下载完成后先别急着分析，做三件最基础的检查（每件都是 R 一行命令）：

```r
# 1. 表达矩阵维度：行数（基因/探针）× 列数（样本）
dim(expr)

# 2. 缺失值数量：缺失太多说明数据质量差，需要处理
sum(is.na(expr))

# 3. 样本分组：看看有哪些组别、每组多少样本
#    注意：先 colnames(pdata) 确认分组列名，再替换下面的列名
table(pdata$characteristics_ch1)
```

这三个检查回答三个问题：**数据有多大？缺失多不多？样本怎么分组？** 如果缺失值很多（例如数以万计），第 07 章会教你处理（过滤 / 补缺 / 剔除），这里先记录观察即可。把检查结果写进 `results/` 下的记录文件，作为分析日志的一部分：

```r
# 记录质控结果（示例），供日后核对
writeLines(c(
  paste("探针数:", dim(expr)[1], "样本数:", dim(expr)[2]),
  paste("缺失值数:", sum(is.na(expr)))
), "results/00_QC_note.txt")
```

**两个容易踩的坑**：一是表达矩阵的列名（GSM 编号）与 `pData` 的行名必须一一对应，分组前最好确认顺序一致（`all(colnames(expr) == rownames(pdata))` 返回 TRUE 才放心）；二是分组列里可能有大小写或写法差异（如 tumor 与 Tumor），用 `unique()` 查看所有取值后再写分组逻辑。

另外，质控的阈值没有统一标准，关键是如实记录并解释你的选择——这本身就是严谨的做法。

## 6.5 数据登记意识

做科研要能“说清楚数据从哪来”。建议在 `data/README.md`（或 data 目录下一个文本文件）里做登记：

```markdown
# 数据登记表
- TCGA-LIHC：来源 UCSC Xena；表达：HTSeq-FPKM；临床：phenotype；
  accession/版本：TCGA-LIHC（以 Xena 页面标注为准）；下载日期：2025-XX-XX
- GSE14520：来源 GEO；平台：GPL3921；
  accession：GSE14520；下载日期：2025-XX-XX
```

只要记住三要素——**accession 编号、下载日期、版本/平台**——以后写论文方法学、或被审稿人追问数据来源时，都能答得上来。（一句话登记就够了，但必须记。）

**为什么登记很重要？** 公共数据库的版本会更新（例如 GEO 偶尔补充样本、Xena 更新注释版本），如果不记录下载日期和版本，以后重跑分析可能得到与论文不一致的结果。另外，审稿人或读者要求“提供数据来源”时，登记表就是你的证据。

## 本章小结

- TCGA-LIHC 用 UCSC Xena 网页下载（HTSeq-FPKM 表达 + phenotype 临床）；GDC 是更权威但更繁琐的备选；
- GSE14520 用 `getGEO("GSE14520", GSEMatrix = TRUE)` 一行下载，`exprs()` 取表达矩阵、`pData()` 取临床信息；探针 ID 需用平台注释（GPL3921）映射为基因 Symbol；
- 原始数据放 `data/raw/`（只读），脚本、结果、图分开存放；
- 下载后先做三维度检查：`dim()`、`sum(is.na())`、`table(分组列)`；
- 登记 accession、下载日期、版本，是数据管理的底线。

## 思考题

1. 为什么表达矩阵的列名是 GSM 编号而不是患者编号？患者编号信息一般藏在 `pData` 的哪一列？
2. 如果 `getGEO("GSE14520")` 返回的 `gse` 有多个元素，说明什么？用 `gse[[2]]` 取出来看看与第一个有什么不同。
3. “探针 ID” 和 “基因 Symbol” 有什么区别？为什么分析前必须做映射？
4. 在 `pData(gse[[1]])` 中找一找 GSE14520 的生存相关字段（如生存时间、生存状态），想一想后续 KM 生存分析需要哪两列。
