从FPKM到Counts：手把手教你准备DESeq2所需的输入数据（附格式转换脚本）
在RNA-seq数据分析中，DESeq2因其稳健的统计模型和出色的差异表达分析能力，已成为生物信息学领域的黄金标准工具。然而，许多初学者在兴奋地安装好DESeq2包后，往往会遇到一个令人沮丧的"入门墙"——工具要求输入原始计数数据（raw counts），但手头只有FPKM或TPM格式的表达矩阵。这种标准化后的数据虽然便于样本间比较，却丢失了DESeq2建模所需的关键统计特性。
本文将系统解决这个"最后一公里"问题，带您完成从常见标准化格式到原始计数的完整转换流程。不同于简单的代码堆砌，我们会深入探讨不同转换方法的数学原理和适用场景，并提供可复用的R脚本和命令行工具，让您无论从GEO数据库下载数据还是处理自有测序结果，都能快速获得DESeq2-ready的输入文件。
1. 为什么DESeq2坚持使用原始计数？
计数数据的统计优势
DESeq2的负二项分布模型直接依赖于原始读段计数的离散特性。当您使用FPKM/TPM时，三个关键信息已经丢失：
- 基因长度偏差（已通过外显子长度标准化消除）
- 文库大小差异（已通过总量标准化消除）
- 计数数据的离散-均值关系（方差随均值变化的特性）
注意：TPM虽然改进了FPKM的跨样本可比性，但依然不适合直接作为DESeq2输入
标准化格式的潜在风险
下表对比了不同数据格式对差异分析的影响：
| 数据格式 | 适合DESeq2 | 保留离散特性 | 保持文库大小信息 | 适用场景 |
|---|---|---|---|---|
| Raw counts | ✓ | ✓ | ✓ | 差异表达分析 |
| FPKM | ✗ | ✗ | ✗ | 样本间表达比较 |
| TPM | ✗ | ✗ | ✗ | 样本间表达比较 |
| RSEM输出 | 有条件✓ | 部分保留 | ✓ | 当包含expected counts时 |
2. 数据溯源：获取原始计数的四种途径
2.1 从公共数据库直接下载count矩阵
GEO/SRA等数据库的提交者有时会同时提供原始计数数据，查找技巧：
- 在数据集页面搜索"count matrix"、"raw counts"等关键词
- 检查Supplementary Files中的
*counts*.txt或*reads*.csv文件 - 对于GSE编号，尝试在GEO查询框输入：
GSEXXXX AND "count"
2.2 使用salmon/kallisto的定量结果
现代准比对工具输出的转录本水平数据可通过tximport转换为基因水平计数：
library(tximport) files <- c("sample1.quant", "sample2.quant") # salmon输出目录 tx2gene <- read.csv("tx2gene.csv") # 转录本到基因的映射表 txi <- tximport(files, type="salmon", tx2gene=tx2gene) dds <- DESeqDataSetFromTximport(txi, ...)2.3 从BAM文件重新计数
当只有比对文件时，使用featureCounts工具生成计数矩阵：
featureCounts -a annotation.gtf -o counts.txt -T 8 -t exon -g gene_id *.bam2.4 FPKM/TPM回算原始计数（最后选择）
当前三种方法都不可行时，可通过以下公式反向估算：
counts ≈ (FPKM × gene_length × total_reads) / 10^9对应的R函数实现：
fpkm_to_counts <- function(fpkm, gene_length, total_reads=1e6) { counts <- fpkm * gene_length * total_reads / 1e9 round(counts) }提示：回算法存在较大误差，仅建议作为最后手段使用
3. 实战演练：GSE数据集处理全流程
以GSE123456数据集为例，演示完整处理流程：
3.1 数据下载与格式识别
library(GEOquery) gse <- getGEO("GSE123456") # 获取元数据 supp_files <- getGEOSuppFiles("GSE123456") # 下载补充文件常见的FPKM文件识别特征：
- 列名包含"FPKM"或"TPM"
- 数值范围为0到数十万
- 行名为基因Symbol或Ensembl ID
3.2 使用tximport处理kallisto输出
当数据集提供SRA原始数据时：
library(readr) samples <- read_csv("samples.csv") # 包含SRR编号和分组信息 files <- file.path("kallisto_output", samples$run, "abundance.h5") names(files) <- samples$run txi <- tximport(files, type="kallisto", txOut=TRUE)3.3 计数矩阵的完整性检查
合格的DESeq2输入矩阵应满足：
- 全部为整数
- 无缺失值
- 行名唯一
- 列名与样本信息表一致
验证代码片段：
stopifnot(all(round(mtx) == mtx)) # 检查整数 stopifnot(!any(is.na(mtx))) # 检查缺失值 stopifnot(nrow(mtx) == length(unique(rownames(mtx)))) # 检查行名唯一性4. 高级技巧与疑难排解
4.1 处理基因ID不一致问题
当行名格式不匹配时，使用biomaRt进行转换：
library(biomaRt) ensembl <- useMart("ensembl", dataset="hsapiens_gene_ensembl") ids <- getBM(attributes=c("ensembl_gene_id", "hgnc_symbol"), filters="hgnc_symbol", values=rownames(mtx), mart=ensembl)4.2 多批次数据的合并处理
合并不同来源的计数矩阵时需注意：
- 确保使用相同的基因注释版本
- 检查批次效应（使用PCA初步判断）
- 考虑使用combat或sva包校正批次
4.3 内存优化策略
对于大型矩阵：
library(HDF5Array) seed <- writeHDF5Array(mtx, "counts.h5") # 将矩阵存储在磁盘 dds <- DESeqDataSetFromMatrix(seed, ...) # 从磁盘懒加载5. 完整转换脚本示例
以下是一个自动化处理FPKM转counts的R脚本框架：
#!/usr/bin/env Rscript # FPKM_to_DESeq2.R - 转换FPKM矩阵为DESeq2输入 suppressPackageStartupMessages({ library(tidyverse) library(optparse) }) option_list <- list( make_option(c("-i", "--input"), help="Input FPKM matrix"), make_option(c("-l", "--length"), help="Gene length file"), make_option(c("-o", "--output"), default="counts.csv", help="Output file") ) args <- parse_args(OptionParser(option_list=option_list)) fpkm2counts <- function(fpkm_file, length_file, output) { # 读取输入文件 fpkm <- read.csv(fpkm_file, row.names=1, check.names=FALSE) gene_len <- read.csv(length_file, row.names=1)[,1] # 检查基因一致性 common_genes <- intersect(rownames(fpkm), names(gene_len)) fpkm <- fpkm[common_genes,] gene_len <- gene_len[common_genes] # 计算总读长（假设中位数FPKM=1对应30M reads） lib_size <- median(colSums(fpkm)) * 1e6 / 30 # 转换计算 counts <- round(fpkm * gene_len * lib_size / 1e9) # 输出结果 write.csv(counts, output) message("Saved counts matrix to ", output) } fpkm2counts(args$input, args$length, args$output)使用方式：
Rscript FPKM_to_DESeq2.R -i fpkm_matrix.csv -l gene_lengths.csv -o deseq2_counts.csv在实际项目中，我们通常会遇到各种特殊情况和数据异常。比如最近处理的一个肿瘤数据集就出现了FPKM值全为零的基因，后来发现是注释版本不匹配导致的。这种情况下，直接剔除这些基因比强行转换更有意义——毕竟DESeq2也会过滤低表达基因。
