RNA-seq下游TPM、FPKM分析方法及ensembl ID-symbol的转化和featureCounts合并

或许各位生信人接过这样的任务——师兄让你把原始数据跑完上游并分析出TPM、FPKM结果交给他？这些内容是简单的下游分析，但是实操起来可能遇到很多问题，又很难找到相关的经验贴，这里旨在给各位提供一套简单可行的流程。首先介绍多批次跑出的count文件如何进行拼接分析。

featureCounts合并方法

首先一定要确保上游分析用了相同的参考基因组和注释文件，合并需要用到R，脚本如下

# 将本次涉及的样本名按文件及列名顺序写下
sample_names <- c("FPpePEP1", "FPpePEP2", "FPpePEP3", "FPpePEP4",
"FS4D51", "FS4D52", "FS4D53",
"H9S4D51", "H9S4D52",
"HPP15", "FPP_P3", "FPP_10D_P2", "FPP_10D_P0_531", "FPP_10D_P0_207")

# 获取文件列表
files <- c("featureCounts_result_1.txt", "featureCounts_result_2.txt")

# 循环读取，为了防止表头出错跳过了说明行和表头行，用上面手写的样品名
count_list <- lapply(files, function(f) {
df <- read.table(f, header=FALSE, skip=2, row.names=1)
return(df[, 6:ncol(df)])
})

# 横向合并所有矩阵
final_matrix <- do.call(cbind, count_list)

# 绑定列名
colnames(final_matrix) <- sample_names

# 1. 检查维度：应该是 (基因数) x 14
dim(final_matrix)

# 2. 打印头部：确保 ENSG ID 在最左侧（行名），样本名在最上方（列名）
head(final_matrix)

# 保存为制表符分隔的文本文件
write.table(final_matrix, file = "combined_counts_final.txt",
sep = "\\t", quote = FALSE, row.names = TRUE, col.names = TRUE)

# 为了方便 EXCEL 打开，也可以保存为 CSV
write.csv(final_matrix, file = "combined_counts_final.csv",
quote = FALSE, row.names = TRUE)

用Gemini写了个自动读取拆分表头的，R基础比较好的朋友可以看看，我样本少，选择用上面那个从count1抄到count2从左往右抄

# ==============================================================================
# CSDN 备选方案：全自动多样本 Count 矩阵合并与列名（长路径）清洗脚本
# 适用场景：样本量极多、不想手动抄写样本名、且 featureCounts 输出带路径时
# ==============================================================================

# 1. 定义需要合并的 featureCounts 结果文件列表
files <- c("featureCounts_result_1.txt", "featureCounts_result_2.txt")

# 2. 循环读取并提取 Count 列
count_list <- lapply(files, function(f) {
# 【关键点 1】header = TRUE：R 会自动跳过开头的 '#' 注释行，并把第二行识别为表头
# 【关键点 2】check.names = FALSE：极其重要！防止 R 把路径中的 '/' 和 '-' 自动转换成 '.'
df <- read.table(f, header = TRUE, row.names = 1, check.names = FALSE)

# 提取表达量列：由于 row.names=1 移走了第一列，原本第 7 列的第一个样本现在平移到了第 6 列
return(df[, 6:ncol(df), drop = FALSE])
})

# 3. 横向合并所有样本矩阵
final_matrix <- do.call(cbind, count_list)

# 4. 【核心自动化】一键提取并清洗长路径列名
# 提取合并后的原始列名（此时带有长路径，如 "/mnt/d/wzc_rnaseq/FPpePEP1.bam"）
raw_names <- colnames(final_matrix)

# 步骤 A：利用 basename() 自动砍掉前面的长路径，只留下文件名 "FPpePEP1.bam"
short_names <- basename(raw_names)

# 步骤 B：利用 gsub() 去掉文件后缀（根据实际情况修改，如 ".bam" 或 ".sorted.bam"）
clean_names <- gsub(".bam", "", short_names)

# 将清洗干净、整洁的样本名重新赋给矩阵列名
colnames(final_matrix) <- clean_names

将脚本置于Count文件相同文件夹，使用R脚本的方法有两种：

# 在Linux系统内
Rscript xxx.R

# 在R内部
source("xxx.R")

Ensembl ID到gene name的转化

由于count、TPM、FPKM文件常进行基于经验的人工分析和筛选，我们需要通俗的基因名便于查询，而不希望看到满文档形如ENSG00000142611的ID。如果大家通过biomaRt 或 org.Hs.eg.db 处理数据的话可能碰到兼容性问题、数据库网络问题和多种报错，这里为大家介绍手动下载的方法，这样下载的txt文件可以永久保存调用。首先我们打开https://www.ensembl.org/数据库，点击上方 "BioMart" 工具；在"-choose dataset-"处选择"Ensembl Genes 116"（可能会更新，总之选Genes）；我的例子是人类，选择Human Genes（GRCh38.p14），如图1；左边需要选的是Attributes，点"Genes"前面的加号；只选"Gene stable ID"和"Gene name"；点做上方"Results"如图2的画面再选右上方"go"可以存为 "mapping.txt" 便于按我后续流程走。

这里参考了Spirallock大佬的文章，后面会注明，我希望这种教程可以被更多人看到

Length的获取

由于后续计算需要用到L即基因的碱基数，我们需要从Count中同时提取这个数据，为了防止出错同时方便后续的阅读，我喜欢另存为gene_lengths.txt文件，脚本如下，比较简单

# 设置文件，有多个 count 的情况下用一个文件就行
file <- c("featureCounts_result_1.txt")

# 这一步只做一次，用于生成 gene_lengths.txt
first_file <- read.table(file, header=FALSE, skip=2, row.names=1)
# 提取行名（GeneID）和 Length（第 6 列）
gene_lengths <- first_file[, 5, drop=FALSE] # 这是因为有一列成 name 了
colnames(gene_lengths) <- c("Length")

# 保存 Length 文件 (用于后续计算 TPM/FPKM)
write.table(gene_lengths, "gene_lengths.txt", sep="\\t", quote=FALSE, row.names=TRUE)

后面有空我也可以补充讲解一下TPM、FPKM的算法

R实现TPM、FPKM输出

做组学的朋友常用R，因此使用简单的R语言脚本实现任务，同时输出Counts、TPM、FPKM文件

# 1. 读取数据
counts <- read.table("combined_counts_final.txt", header=TRUE, row.names=1)
lengths <- read.table("gene_lengths.txt", header=TRUE, row.names=1)
map <- read.table("mapping.txt", header=TRUE, sep="\\t", stringsAsFactors=FALSE)

# 2. 计算 TPM 和 FPKM
rpk <- counts / (lengths$Length / 1000)
tpm <- t(t(rpk) / colSums(rpk) * 1e6)
fpkm <- t(t(counts / (lengths$Length / 1000)) / (colSums(counts) / 1e6))

# 3. 基因 ID 映射与唯一性处理
# 将 map 的第一列（ID）设为行名，方便直接用矩阵索引提取
rownames(map) <- map[, 1]

# 根据表达矩阵的行名，提取对应的 Gene name（第二列）
all_ensg <- rownames(counts)
gene_names <- map[all_ensg, 2]

# 逻辑替换：如果没匹配到 Symbol（为 NA 或空字符串），则填补回原始的 Ensembl ID
gene_names <- ifelse(is.na(gene_names) | gene_names == "", all_ensg, gene_names)

# 去重逻辑：利用 make.unique 自动给重复的 Symbol 加上后缀（如 .1, .2），确保矩阵行名唯一
unique_gene_names <- make.unique(as.character(gene_names))

# 4. 保存为 CSV 文件的函数
save_as_csv <- function(mat, names, filename) {
df <- as.data.frame(mat)
rownames(df) <- names
write.csv(df, filename, quote=FALSE)
}

# 执行保存：同时导出 TPM、FPKM 和原始 Counts 的 Symbol 版本
save_as_csv(tpm, unique_gene_names, "TPM_Symbol.csv")
save_as_csv(fpkm, unique_gene_names, "FPKM_Symbol.csv")
save_as_csv(counts, unique_gene_names, "Counts_Symbol.csv")

message("处理完成！已自动处理重复名称，成功生成 TPM、FPKM 和 Counts 的 CSV 文件。")

参考文献：

python实现Ensembl ID和gene symbol的相互转换