肝癌耐药机制探索：基于deseq2、edgeR与limma_voom的转录组差异分析实战

原创

于 2026-03-06 00:46:46 发布 · 352 阅读

1. 从肝癌耐药说起：为什么转录组差异分析是关键第一步

大家好，我是老张，在生物信息这个圈子里摸爬滚打了十几年，主要跟各种组学数据和算法打交道。今天想和大家深入聊聊一个在肿瘤研究，特别是耐药机制探索中，几乎人人都会遇到，但又常常让人头大的技术环节——转录组差异分析。

就拿肝癌来说，索拉菲尼这类靶向药一开始效果可能不错，但用着用着，癌细胞就“学聪明了”，产生了耐药性。这背后的分子机制是什么？是哪些基因在“捣鬼”？它们的表达量是升高了还是降低了？要回答这些问题，我们第一步要做的，就是从海量的转录组测序数据里，把那些在耐药细胞和敏感细胞之间表达有显著差异的基因给“揪”出来。这个过程，就是差异表达分析。

你可能会问，现在做差异分析的R包那么多，DESeq2、edgeR、limma-voom，我该用哪个？网上教程一搜一大把，但照着代码跑一遍，结果出来了，心里却更没底了：这几个方法得出的结果怎么不太一样？哪个更可靠？最后发文章到底该用哪个列表的基因？这些问题，我当年刚开始做分析的时候也纠结了很久，踩过不少坑。

所以，这篇文章我不打算只给你扔几段代码。我想结合一个真实的肝癌耐药数据集（GSE213615），带你完整走一遍从数据下载、预处理，到用三种主流方法（DESeq2, edgeR, limma-voom）分别做分析，再到最后如何交叉验证、解读结果的实战流程。我的目标是，你读完不仅能自己跑通流程，更能理解每个步骤背后的“为什么”，知道在不同情况下如何选择和判断，真正把工具用活，为你的耐药机制研究打下坚实可靠的第一步。

2. 实战起手式：数据获取与预处理，坑都在这儿了

工欲善其事，必先利其器。分析的第一步，是把数据干干净净地准备好。这一步看似琐碎，却至关重要，很多后续分析出的怪问题，根源都在这儿。我们用的数据集是GEO数据库里的GSE213615，这个数据研究的是两种肝癌细胞系（HepG2和Huh7）经索拉菲尼处理后产生的耐药细胞，对比未处理的对照组，非常适合我们今天的主题。

2.1 数据下载与初步整理

首先，我们需要从GEO下载原始数据。这里我强烈推荐使用GEOquery这个R包，它是和GEO数据库交互的“瑞士军刀”。原始文章里提供的代码是从本地已下载的文件读取，为了更通用，我先演示如何直接从GEO获取。

# 加载必要的包
library(GEOquery)
library(dplyr)

# 指定GSE编号
proj <- "GSE213615"

# 下载数据集，这可能会花点时间，取决于网速
gset <- getGEO(GEO = proj, GSEMatrix = TRUE, getGPL = FALSE, destdir = ".")

# 通常getGEO会返回一个列表，我们取第一个元素
gset <- gset[[1]]

# 提取表达矩阵（这个数据集作者已经上传了处理过的count矩阵）
expr <- exprs(gset)
# 提取样本临床信息（表型数据）
clinical <- pData(gset)

# 先看一眼数据长什么样
dim(expr) # 查看基因数和样本数
expr[1:5, 1:3] # 看看前5个基因在前3个样本的表达值
head(clinical[, 1:5], 3) # 看看临床信息的前几列

下载后你会发现，这个数据集的表达矩阵可能不是原始的“count”数据，而是经过了一些标准化处理。这里是一个超级重要的点：DESeq2和edgeR这两个包是专门为基于计数的转录组数据设计的，它们要求输入是原始的整数read counts。 如果输入是FPKM、TPM等标准化后的连续值，结果会不准确。好在GSE213615提供了每个样本单独的原始count文件，我们需要像原始文章那样，把它们合并起来。

2.2 合并原始Count文件与基因过滤

假设你已经把从GEO Supplementary里下载的多个样本的count文件放到了同一个文件夹下。接下来的操作就像是一个数据清洗流水线：

# 设置存放原始count文件的目录
file_directory <- "./GSE213615_RAW"

# 列出目录下所有文件
fs <- list.files(file_directory, pattern = "*.txt|*.gz", full.names = FALSE)

# 使用循环或lapply读取并合并，这里提供一个更易理解的循环版本
exp_list <- list() # 创建一个空列表存放每个样本的数据

for (i in 1:length(fs)) {
  file_path <- file.path(file_directory, fs[i])
  # 读取文件，注意分隔符可能是制表符\t
  dat <- read.table(file_path, header = TRUE, sep = "\t", row.names = NULL, stringsAsFactors = FALSE)

  # 数据清洗：去除描述列中含有“lncRNA”的行（根据原始数据说明）
  dat <- dat[!grepl("lncRNA", dat$description), ]

  # 提取我们需要的列：基因Symbol列，以及样本表达量列（列名包含Hep或Huh）
  symbol_col <- which(colnames(dat) == "symbol")
  sample_cols <- which(grepl("Hep|Huh", colnames(dat))) # 匹配列名

  # 如果找不到样本列，可能需要调整匹配规则，这里只是示例
  if(length(sample_cols) > 0){
    dat <- dat[, c(symbol_col, sample_cols)]
  } else {
    # 如果列名不匹配，可能需要查看原始列名
    print(head(colnames(dat)))
    stop("未找到包含'Hep'或'Huh'的样本列，请检查数据。")
  }

  # 从文件名中提取样本ID（例如GSM6589876）
  sample_id <- strsplit(fs[i], "_")[[1]][1]
  # 将唯一的样本ID设为列名（除了symbol列）
  colnames(dat)[-1] <- sample_id # 假设第一列是symbol

  # 存入列表
  exp_list[[sample_id]] <- dat
}

# 使用reduce或循环，基于symbol列合并所有数据框
# 这里使用dplyr的full_join进行迭代合并
library(purrr)
exp_merged <- reduce(exp_list, full_join, by = "symbol")

# 处理缺失值：对于合并后某些样本缺失的基因，计数填充为0（表示未检测到）
exp_merged[is.na(exp_merged)] <- 0

# 将gene symbol设为行名
rownames(exp_merged) <- exp_merged$symbol
exp_merged$symbol <- NULL # 移除多余的symbol列

# 查看合并后的数据
dim(exp_merged)
exp_merged[1:5, 1:5]

数据合并好了，但还不能直接用。转录组数据里有很多基因在所有或大多数样本里表达量极低或为零，这些基因没有统计分析的效力，留着只会增加计算负担和多重检验校正的压力。所以我们需要做基因过滤。

# 过滤低表达基因：保留在至少一半样本中表达量大于0的基因
# 注意：这里的“表达量大于0”是针对count数据。如果是非常深度的测序，有时会用大于1或5作为阈值。
keep_genes <- rowSums(exp_merged > 0) >= 0.5 * ncol(exp_merged)
table(keep_genes) # 看看过滤掉了多少基因
exp_filtered <- exp_merged[keep_genes, ]

# 确保数据是整数（count数据的本质）
exp_filtered <- round(exp_filtered)
dim(exp_filtered)