2. 数据获取与质量控制

好，咱们开始进入实操环节。做多组学分析，第一步就是找数据。你想想看，没有高质量的数据，后面那些花里胡哨的模型全是白搭。所以这一节，我重点聊聊三个最常用的公开数据库——TCGA、GEO、ENCODE，以及怎么把数据弄下来、怎么判断数据能不能用。

2.1 三大公开数据库，你总得用熟一个

做精准医学研究，离不开这几个数据宝库。我个人的习惯是：先看TCGA有没有你需要的癌种，没有就去GEO淘，需要调控元数据就找ENCODE。咱们一个个说。

2.1.1 TCGA（癌症基因组图谱）

TCGA是迄今为止最系统的癌症多组学数据库。它收集了33种癌症类型，超过2万个样本的基因组、转录组、表观遗传组、蛋白质组数据。说白了，你想研究哪个癌种的突变和表达，TCGA基本都能覆盖。

核心数据类别包括：

• 体细胞突变（MAF格式）
• mRNA表达（RNA-Seq，FPKM/TPM）
• miRNA表达
• DNA甲基化（450K/27K芯片）
• 拷贝数变异（CNV，GISTIC处理结果）
• 临床数据（生存信息、分期、治疗史）

2.1.2 GEO（基因表达综合数据库）

GEO的数据就杂多了。它不光包含癌症，还有正常组织、发育、疾病状态、药物处理等各种条件。我遇到过最头疼的问题就是——数据格式不统一。GEO里有芯片数据、有测序数据、甚至有甲基化芯片，整理起来特别花时间。

GEO的核心结构是：

层级	含义	例子
GSE系列	一个完整的研究项目	GSE12345
GSM样本	单个样本的实验结果	GSM123456
GPL平台	检测使用的芯片或测序平台	GPL570 (Affymetrix)

2.1.3 ENCODE（DNA元件百科全书）

ENCODE的核心任务是注释人类基因组中的功能元件。如果TCGA和GEO告诉你基因表达什么样，ENCODE就告诉你为什么这样——通过转录因子结合位点、组蛋白修饰、染色质可及性（ATAC-seq）等数据。

常用数据类型：

• ChIP-seq（转录因子、组蛋白修饰标记）
• DNase-seq / ATAC-seq（染色质开放区域）
• RNA-seq（不同细胞系表达谱）
• Hi-C（三维基因组结构）

嗯，做多组学整合分析时，ENCODE的数据常常用来构建调控网络。比如你想分析某个转录因子在TCGA的乳腺癌样本中是否异常活跃，就可以用ENCODE的ChIP-seq峰来定义TF的靶基因集。

2.2 数据下载，别再用浏览器手动点了

很多初学者还在网页上一个个点下载，其实有专门的工具。我自己的工作流基本不离开命令行。

2.2.1 TCGA数据下载：gdc-client

TCGA现在全归到GDC（Genomic Data Commons）管理了。官方提供了命令行工具gdc-client，支持断点续传和批量下载。

常用命令：

# 下载manifest文件里列出的所有文件
gdc-client download -m gdc_manifest.txt

# 按项目筛选下载（比如只下载BRCA）
gdc-client download -m gdc_manifest.txt --dir ./BRCA_data

这里有个坑：gdc-client默认下载的是原始BAM或VCF文件，大小动辄几十GB。如果你只想做差异表达分析，建议下载RNA-Seq的STAR计数结果（Gene Expression Quantification），文件小很多。

2.2.2 GEO数据下载：GEOquery + wget

GEO的情况稍微复杂点。我一般用R包GEOquery来拉元数据和表达矩阵。对于高通量测序数据，GEO通常在SRA里存放原始序列，需要用fastq-dump或prefetch来获取。

# R语言示例
library(GEOquery)

# 下载系列矩阵文件
gse <- getGEO("GSE12345", destdir = ".")

# 提取表达矩阵
exprs_data <- exprs(gse[[1]])

# 下载原始测序数据（需要安装SRA Toolkit）
prefetch SRR1234567
fastq-dump --split-files SRR1234567

你想想看，如果数据量不大（比如几十个样本的芯片数据），直接在网页下载MINiML文件也够用了。但要是几百个RNA-Seq样本，我还是建议用prefetch配合SRA。

2.2.3 ENCODE数据下载：API调用

ENCODE提供了RESTful API，可以编程式批量获取metadata和文件链接。我用Python写过一个小脚本，自动根据实验类型、生物样本、文件格式来筛选下载。

# Python示例：通过ENCODE API搜索并下载ChIP-seq峰
import requests
import json

url = "https://www.encodeproject.org/search/?type=Experiment"
params = {
    "assay_title": "ChIP-seq",
    "biosample_ontology.classification": "cell line",
    "target.label": "CTCF",
    "status": "released",
    "format": "json"
}
resp = requests.get(url, params=params)
data = resp.json()

# 获取第一个实验的bed文件链接
for graph in data["@graph"][0]["files"]:
    if graph["file_type"] == "bed narrowPeak":
        download_url = graph["s3_uri"]
        print(f"Download: {download_url}")
        break

这个API返回的数据结构很规范，我建议多看几次返回的JSON，搞清楚@graph、files、replicates这些字段含义，写起代码来就顺手多了。

2.3 质量评估，不合格的数据不如不用

数据下载下来，千万别直接开跑。我见过最惨的案例是——用了一个测序深度严重不足的样本做差异分析，结果发现差异表达基因全是低表达噪音。质量控制这几步，省不了。

2.3.1 测序数据质量评估：FastQC

对于RNA-Seq、ChIP-Seq等测序数据，第一步是跑FastQC。它会生成每个样本的质量报告，包括：碱基质量分布、GC含量、接头污染、重复比例等。

# 创建输出目录
mkdir -p fastqc_reports

# 运行FastQC
fastqc -o fastqc_reports/ sample1_R1.fastq.gz sample1_R2.fastq.gz

# 批量跑
ls *.fastq.gz | parallel -j 4 "fastqc -o fastqc_reports/ {}"

重点关注几个指标：

• Per base sequence quality：整个读长范围内碱基质量在Q30以上（错误率<0.1%）
• Adapter Content：接头残留比例应接近0%
• Overrepresented sequences：高比例重复序列可能表明样本污染

2.3.2 表达谱质量评估：PCA + 热图

拿到表达矩阵后，我做的第一件事不是差异分析，而是看PCA。目的：检查是否有异常样本、批次效应是否明显。

# R示例：PCA诊断
library(ggplot2)
library(DESeq2)

# 假设count_matrix是基因×样本的计数矩阵，colData包含样本分组信息
dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = colData,
                              design = ~ condition)

# 使用vst或rlog变换后跑PCA
vsd <- vst(dds, blind = TRUE)
pca_data <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)

# 画图看看样本是否按预期分组
ggplot(pca_data, aes(PC1, PC2, color = batch, shape = condition)) +
    geom_point(size = 3) +
    theme_minimal()

如果某个样本远离整个群体，我通常会查看它的测序深度、比对率、基因检出数。要是这些指标也异常，这个样本就得准备扔了。

2.3.3 常见质量控制指标一览

数据类型	核心指标	合格标准	常见问题指标
RNA-Seq	比对率	>70%到参考基因组	比对率低→可能样本降解
RNA-Seq	基因检出数	通常>10000个基因（FPKM>1）	检出数过低→测序深度不够
ChIP-Seq	FRiP	峰内读段比例>1%	FRiP<0.5%→抗体或IP失败
WGS/WES	平均深度	WGS>30x, WES>80x	深度不足→低频突变漏检
DNA甲基化（450K）	检测p值	每个探针p值<0.05	探针未检出→杂交或扩增失败

2.4 本章知识体系总览

下面这张图总结了我讲的数据获取与质量控制流程，从原始数据资源到最终可用的清洗数据，每一步都不能少。