第3章 RNA-seq数据获取：TCGA与GEO数据库下载流程

做肿瘤亚型分类，第一步永远不是调模型，而是搞到对的数据。我见过不少同学拿到公开数据直接用，结果跑完才发现样本批次对不上、临床信息不全、甚至基因ID格式都不统一——那真是浪费时间。

这一章我们就来解决数据源头的问题。说白了，就是要搞清楚两个最大的癌症基因组学数据仓库：TCGA和GEO，它们的RNA-seq数据怎么拿、拿什么、拿到后怎么管。

3.1 两大数据库：TCGA vs GEO

这俩名字你肯定听过，但它们到底有什么区别？我记刚入行时也迷糊过一阵子。直接看表吧：

对比项	TCGA	GEO
全称	The Cancer Genome Atlas	Gene Expression Omnibus
数据规模	~20,000 原发性癌症样本，33种癌型	超过 200,000 个系列，包含各种疾病
数据类型	RNA-Seq, WXS, WGS, 甲基化, 临床	表达芯片, RNA-Seq, 甲基化, ChiP-Seq等
数据标准	统一处理流程（标准化程度高）	各个实验室自行上传，格式满天飞
适用场景	肿瘤亚型分类、预后标志物、多组学	差异表达、小样本验证、非肿瘤研究
下载方式	GDC API, GDC Data Portal, TCGAbiolinks	GEO Query, GEOmetadb, NCBI SRA Toolkit

3.2 TCGA数据下载：两条主流路线

TCGA数据怎么拿？嗯，这里我走过不少弯路。目前最靠谱的有两条路：

路线一：GDC Data Portal + 手动下载（适合小规模）

去 portal.gdc.cancer.gov，点“Repository”。然后按以下步骤筛选：

1. Cases → Program: TCGA → Project: 选你需要的癌型（比如TCGA-BRCA乳腺癌）
2. Files → Data Category: Transcriptome Profiling → Data Type: Gene Expression Quantification
3. Files → Workflow Type: STAR - Counts（这是官方推荐的定量流程）
4. 点“Add All Files”到购物车，然后下载 Manifest 或者直接下载 Cart。

但说实话，手动下载几百个样本搞死人。所以我一般走路线二。

路线二：GDC API + Python批量下载（推荐）

我写过一个小脚本，用 requests 和 gdc-api-client 来拉数据。下面给个核心模板：

import requests
import json
import os

# 1. 查询TCGA-BRCA的RNA-Seq count文件
query_url = "https://api.gdc.cancer.gov/files"
params = {
    "filters": {
        "op": "and",
        "content": [
            {"op": "in", "content": {"field": "cases.project.project_id", "value": ["TCGA-BRCA"]}},
            {"op": "in", "content": {"field": "files.data_category", "value": ["Transcriptome Profiling"]}},
            {"op": "in", "content": {"field": "files.data_type", "value": ["Gene Expression Quantification"]}},
            {"op": "in", "content": {"field": "files.experimental_strategy", "value": ["RNA-Seq"]}}
        ]
    },
    "fields": "file_id,file_name,cases.submitter_id",
    "format": "JSON",
    "size": 1000
}

response = requests.post(query_url, json=params)
file_list = response.json()["data"]["hits"]

# 2. 下载每个文件
data_endpoint = "https://api.gdc.cancer.gov/data"
for f in file_list[:5]:   # 先下载5个测试
    file_id = f["file_id"]
    file_name = f["file_name"]
    response = requests.get(f"{data_endpoint}/{file_id}", stream=True)
    with open(file_name, "wb") as out:
        for chunk in response.iter_content(chunk_size=1024*1024):
            out.write(chunk)
    print(f"Downloaded: {file_name}")

3.3 GEO数据下载：从系列到表达矩阵

GEO 的数据结构比 TCGA 乱多了。一个GEO系列（GSE）下可能有多个平台（GPL），每个平台对应不同的探针。我总结了一套标准流程：

步骤1：找到目标GSE号

去 ncbi.nlm.nih.gov/geo，搜索关键词加癌型，比如 "breast cancer RNA-Seq"。一般优先挑 Expression profiling by high throughput sequencing 的系列。

步骤2：用GEOquery下载表达矩阵

这是R语言派的经典工具，但我喜欢用Python替代——GEOparse 包很好用：

pip install GEOparse

import GEOparse

# 下载GSE系列数据
gse = GEOparse.get_GEO(geo="GSEXXXXX", destdir="./geo_data", silent=True)

# 提取表达矩阵
expression = gse.pivot_samples("VALUE")
print(expression.shape)

# 获取临床信息（phenotype data）
pheno = gse.phenotype_data
print(pheno.columns)

GEOparse 会自动缓存压缩文件，第二次拉就不用重新下载了。嗯，这个很省心。

步骤3：处理原始RNA-Seq数据（如果需要）

如果GEO只提供了原始fastq文件，那就得用SRA Toolkit下载然后自己跑定量。我个人建议先看看Series Supplementary Files里有没有预先算好的表达矩阵，没有才走下面这条路：

# 下载SRA文件
prefetch SRRXXXXXX

# 转换为fastq
fastq-dump --split-files SRRXXXXXX.sra

# 然后用salmon/kallisto定量（后面章节会专门讲）
salmon quant -i salmon_index -l A -1 sample_R1.fastq -2 sample_R2.fastq -o quant_output

3.4 数据整理与ID转换：统一基因符号

拿到表达矩阵后，八成会遇到基因ID不统一的问题。TCGA给的是 ENSEMBL基因ID（如 ENSG000001...），GEO给的可能还是 SYMBOL 或者 ENTREZ ID。统一转换是必须的一步。

我的工具链：用 mygene 库批量查：

import mygene
import pandas as pd

mg = mygene.MyGeneInfo()

# 假设你有一列ENSEMBL ID
ensembl_ids = ["ENSG000001.1", "ENSG000002.2"]   # 注意去除版本号
# 去除版本号，mygene默认处理不带点的id
clean_ids = [eid.split(".")[0] for eid in ensembl_ids]

# 批量查询，获得symbol
result = mg.querymany(clean_ids, scopes="ensembl.gene", fields="symbol", species="human")
gene_map = {item["query"]: item.get("symbol", None) for item in result}

# 替换表达矩阵的index
expr.index = expr.index.map(lambda x: gene_map.get(x.split(".")[0], x))

注意，mygene 有调用频率限制，大量ID（超过1000个）建议分批查询，加个几秒的延迟，不然会返回空数据。

3.5 数据质量快速检查（别跳过这一步）

下载完数据，我建议先做个快速QC，别等跑完模型才发现数据有问题：

• 检查缺失值：每个基因的count是否为0？异常值太多可能代表样本污染。
• 看文库大小：总counts量级应该在1e6~3e7之间，差太多要考虑是否标准化。
• 看基因覆盖度：蛋白编码基因数量是否在2万左右？如果只有几千，可能是测序深度不够。
• 检查性别相关基因：比如XIST在女性样本里高表达，在男性里低表达，可以用这个做样本标签校验。

3.6 本章小结：数据获取流程全景

下面这张图把整个流程串起来了，你可以对照着走一遍：

看到这个流程，你是不是觉得也不是那么复杂？其实说白了就是两步：选对数据库，然后用对的工具下载。后面几章我们会拿这些干净的数据做差异表达、特征选择，再到分类模型，一步步走下去。