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处理结果)
- 临床数据(生存信息、分期、治疗史)
小提示:我建议你优先下载经Firehose流程标准化后的数据,比如"stddata"。这些数据做了批次校正,不同癌种之间可比性更好。
2.1.2 GEO(基因表达综合数据库)
GEO的数据就杂多了。它不光包含癌症,还有正常组织、发育、疾病状态、药物处理等各种条件。我遇到过最头疼的问题就是——数据格式不统一。GEO里有芯片数据、有测序数据、甚至有甲基化芯片,整理起来特别花时间。
GEO的核心结构是:
| 层级 | 含义 | 例子 |
|---|---|---|
| GSE系列 | 一个完整的研究项目 | GSE12345 |
| GSM样本 | 单个样本的实验结果 | GSM123456 |
| GPL平台 | 检测使用的芯片或测序平台 | GPL570 (Affymetrix) |
注意:同一GSE下的样本可能使用了不同GPL平台,这时候做差异表达分析必须去除批次效应。我曾经踩过这个坑——直接合并两个平台的表达矩阵,结果发现差异基因全是平台本身的信号。
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),文件小很多。
我的推荐做法:在GDC Data Portal上先用Filters(癌症类型、数据类型、文件格式)筛选好,导出Manifest文件,再用gdc-client批量下载。这样一次搞定,不容易乱。
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:高比例重复序列可能表明样本污染
注意:FastQC报告里的“FAIL”不一定全要排除。比如RNA-Seq样品的重复序列比例高是正常的(核糖体RNA富集)。你需要结合实验设计和数据类型来判断哪个FAIL是致命的。
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 本章知识体系总览
下面这张图总结了我讲的数据获取与质量控制流程,从原始数据资源到最终可用的清洗数据,每一步都不能少。
核心要点:数据获取不是简单的下载,质量控制也不是走过场。我自己的经验是:花在质控上的时间,能省下后面分析时80%的排查精力。每多花10分钟检查一个样本,就可能避免一整个项目的无效重复劳动。