从FPKM到Counts：手把手教你准备DESeq2所需的输入数据（附格式转换脚本）

从FPKM到Counts：手把手教你准备DESeq2所需的输入数据（附格式转换脚本）

在RNA-seq数据分析中，DESeq2因其稳健的统计模型和出色的差异表达分析能力，已成为生物信息学领域的黄金标准工具。然而，许多初学者在兴奋地安装好DESeq2包后，往往会遇到一个令人沮丧的"入门墙"——工具要求输入原始计数数据（raw counts），但手头只有FPKM或TPM格式的表达矩阵。这种标准化后的数据虽然便于样本间比较，却丢失了DESeq2建模所需的关键统计特性。

本文将系统解决这个"最后一公里"问题，带您完成从常见标准化格式到原始计数的完整转换流程。不同于简单的代码堆砌，我们会深入探讨不同转换方法的数学原理和适用场景，并提供可复用的R脚本和命令行工具，让您无论从GEO数据库下载数据还是处理自有测序结果，都能快速获得DESeq2-ready的输入文件。

1. 为什么DESeq2坚持使用原始计数？

计数数据的统计优势
DESeq2的负二项分布模型直接依赖于原始读段计数的离散特性。当您使用FPKM/TPM时，三个关键信息已经丢失：

• 基因长度偏差（已通过外显子长度标准化消除）
• 文库大小差异（已通过总量标准化消除）
• 计数数据的离散-均值关系（方差随均值变化的特性）

注意：TPM虽然改进了FPKM的跨样本可比性，但依然不适合直接作为DESeq2输入

标准化格式的潜在风险
下表对比了不同数据格式对差异分析的影响：

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 *.bam

2.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也会过滤低表达基因。