2、常见质控指标：测序深度、比对率、基因检出数、文库复杂度、批次效应
多组学数据拿到手，第一件事不是跑差异分析，而是先做质控。我见过太多人一上来就怼着p值调，结果跑出来的结果全是批次噪音——说白了，数据质量不行，后面一切白搭。
这一章，我把五个最核心的质控指标拆开揉碎了讲。每个指标我都会告诉你：它到底在衡量什么、什么算好什么算差、以及我在真实项目里踩过的坑。
2.1 测序深度
测序深度，说白了就是你对基因组上某个位点做了多少次独立的测序。深度越高，检测到低表达基因或低频突变的能力就越强。
我个人习惯先问一句：你这个实验到底想抓什么信号？如果是检测稀有克隆，深度不够就是致命伤。
| 应用场景 | 推荐测序深度 | 说明 |
|---|---|---|
| 全基因组测序 (WGS) | 30-60x | 肿瘤样本建议60x以上 |
| 全外显子测序 (WES) | 100-200x | 检测低频体细胞突变需要高深度 |
| RNA-seq | 20-30M reads/样本 | 关注的是表达量而非基因组全覆盖 |
| 单细胞RNA-seq | 每细胞 50k-200k reads | 深度过高可能导致UMI饱和 |
检查深度最简单的方式就是看BAM文件的覆盖深度分布。嗯，这里我常用的命令：
samtools depth -b target_regions.bed sample.bam | \
awk '{sum+=$3; count++} END {print "平均深度:", sum/count}'还可以用mosdepth跑更全面的统计：
mosdepth --by target.bed sample sample.bam
# 输出 sample.regions.bed.gz
zcat sample.regions.bed.gz | \
awk '$4 >= 20 {covered++; total++} END {print "20x覆盖率:", covered/total*100, "%"}'2.2 比对率
比对率，测序数据能成功比对到参考基因组上的read比例。比对率低，一种可能是样本污染，另一种可能是物种不匹配。
我做过一个跨物种项目，当时把人的样本比到了小鼠基因组上，比对率只有30%。你想想看，那时候我还傻乎乎地以为是小鼠样本特殊，查了半天才发现是样本标签贴错了。多组学数据，样本来源确认永远是第一步。
比对率通常用STAR、Hisat2或BWA-MEM跑完后，从日志里直接提取：
# STAR 日志
grep "Uniquely mapped reads %" Log.final.out
# 输出示例:
# Uniquely mapped reads % | 85.42%
# BWA-MEM
samtools flagstat aligned.bam | grep "mapped"
# 输出示例:
# 98765432 + 0 mapped (92.1% : N/A)| 测序类型 | 合格比对率 | 常见低比对原因 |
|---|---|---|
| 人全基因组 | > 95% | 接头污染、低质量碱基 |
| RNA-seq (人类) | > 80% | 核糖体RNA残留、比对到其它物种 |
| 宏基因组 | 20-60% | 宿主DNA污染、参考基因组不完整 |
为什么会有人类样本比对率低于80%？我遇到过一次：样品RNA提取时掺入了大量细菌RNA，结果建库跑出来一大半reads比到了NCBI的细菌数据库。我的解决办法：先对原始fastq做一圈kraken2分类，看看到底是什么东西混进来了。
kraken2 --db /path/to/kraken2_db --threads 8 \
--report kraken_report.txt \
sample_R1.fastq sample_R2.fastq2.3 基因检出数
基因检出数，就是在你的数据中能检测到的基因数量。对于RNA-seq，这代表转录组的覆盖广度。
单细胞数据更敏感。一个细胞如果只检出几百个基因，基本可以判断是死细胞或空液滴。我一般设的阈值是：<500基因/细胞直接丢掉。
如何计算基因检出数？用featureCounts或htseq-count得到基因表达矩阵，然后统计非零表达基因数：
# R 语言快速统计
counts <- read.table("gene_counts.tsv", header=TRUE, row.names=1)
detected_genes <- rowSums(counts > 0)
summary(detected_genes)
# 输出示例:
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 234 5893 10432 11234 15423 28934| 数据类型 | 预期检出基因数 | 异常低值判断 |
|---|---|---|
| 整组织RNA-seq (人) | 10,000-20,000 | < 8,000 需复查 |
| 单细胞RNA-seq (每个细胞) | 1,000-5,000 | < 500 细胞活性有问题 |
| ATAC-seq (peak数) | 50,000-100,000 | < 30,000 信号太弱 |
2.4 文库复杂度
文库复杂度衡量的是你的文库中独特DNA/RNA分子占所有测序读数的比例。说白了，就是你的文库到底有多「多样」。
复杂度低意味着大量reads是PCR重复。你想想看，如果90%的reads都是同一个片段的不同拷贝，你的数据利用率就极低。
计算文库复杂度常用Preseq，或者直接看samtools标记后的重复率：
# 用 picard 标记重复
picard MarkDuplicates \
I=aligned.bam O=marked.bam M=dup_metrics.txt
# 查看重复率
grep "PERCENT_DUPLICATION" dup_metrics.txt
# 输出示例:
# PERCENT_DUPLICATION 0.32预设期望值如下：
| 类型 | 期望重复率 | 复杂度判断 |
|---|---|---|
| WGS (全基因组) | < 10% | 极低的重复是正常的 |
| ChIP-seq | 30-70% | 重复不可避免，但>70%要警觉 |
| RNA-seq | < 30% | 高表达基因会产生重复，但整体不能太高 |
| 单细胞RNA-seq | < 20% | UMI设计本就是为了消除重复 |
2.5 批次效应
批次效应是多组学数据中最隐蔽、最破坏性、也最容易被忽视的问题。两个条件完全相同的样本，只因为不同天处理、不同试剂批号、不同操作员，产生的数据就可能出现系统性差异。这种差异就是批次效应。
我最早接触批次效应是在整合两个公共数据集时。看起来条件分组很明显，但PCA一跑，明显分成上下两簇——不是按处理分，是按数据集分的。当时我一拍脑袋，差点当成生物学差异去解读。其实，批次效应的检查方法并不复杂。
首先，用PCA或t-SNE可视化样本聚类：
# R 代码示例
pca <- prcomp(t(log2(normalized_counts+1)), scale.=TRUE)
plot(pca$x, col=factor(batch_info), pch=19)
# 如果样本按batch分组而非实验组，说明批次效应显著其次，用limma的removeBatchEffect进行校正：
# 校正批次效应
library(limma)
corrected_counts <- removeBatchEffect(
normalized_counts,
batch = sample_info$batch,
design = model.matrix(~ treatment, data = sample_info)
)| 检查方法 | 预期结果 | 异常信号 |
|---|---|---|
| PCA按batch着色 | 样本混合均匀 | 样本按batch聚成独立簇 |
| 样本相关性热图 | 同组内高相关 | 同批次间相关性远高于同组 |
| 差异基因数量 | 批次间差异基因很少 | 批次间差异基因多于处理组 |
五个指标联合审查
单独看一个指标，不足以判断数据质量。我总结了一个快速检查清单，每次拿到新数据都过一遍：
- 测序深度：全局均值达标？局部有没有深度塌陷的hotspot？
- 比对率：整体在合格线上？有没有异常低比对的样本？
- 基因检出数：分布是否合理？有无双峰？
- 文库复杂度：重复率高不高？是不是由特定区域引起？
- 批次效应：PCA/相关性矩阵是否显示出非生物因素的聚集？