R语言实战：从TCGA官网下载到火山图，手把手搞定肝癌(LIHC)差异表达分析全流程
R语言实战：从TCGA官网下载到火山图，手把手搞定肝癌(LIHC)差异表达分析全流程
在生物信息学研究中，TCGA数据库是癌症基因组分析的重要资源宝库。对于肝癌(LIHC)研究而言，掌握从原始数据获取到差异表达分析的全流程，是每个科研人员必备的核心技能。本文将带你用R语言完成从TCGA数据下载、预处理、差异分析到可视化呈现的完整过程，特别针对初学者可能遇到的各类问题提供解决方案。
1. TCGA数据获取与前期准备
1.1 数据下载流程详解
访问TCGA官方数据门户(https://portal.gdc.cancer.gov/)后，按以下步骤操作：
- 在Cohort Builder中选择Program为TCGA，Project为LIHC
- 进入Repository后，在侧边栏设置筛选条件：
- Experimental Strategy: RNA-Seq
- Data Category: Transcriptome profiling
- Data Type: Gene Expression Quantification
- 将筛选结果全部加入购物车(Cart)
- 下载两个关键文件：
- Sample Sheet（样本信息表）
- Cart文件（数据下载清单）
注意：下载的压缩包通常包含数百个文件，建议准备至少10GB的存储空间
1.2 R环境配置与包安装
在开始分析前，需确保R环境中已安装必要工具包：
# 基础数据处理包 install.packages(c("data.table", "dplyr", "stringr")) # 生物信息学专用包 if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c("edgeR", "limma", "DESeq2")) # 可视化包 install.packages(c("ggplot2", "ggrepel", "ggprism"))常见安装问题解决方案：
- 遇到Bioconductor镜像问题，可尝试：
options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor") - 包依赖冲突时，建议新建干净的R session
2. 数据整理与质量控制
2.1 文件解压与目录结构
推荐建立清晰的文件夹体系：
TCGA-LIHC/ ├── RawMatrix/ # 存放原始下载数据 ├── ProcessedData/ # 存放处理后的矩阵 ├── Results/ # 分析结果 └── Scripts/ # R脚本解压数据的R代码：
# 设置工作路径 setwd("~/TCGA-LIHC") # 解压下载包 untar("gdc_download_20240606_135942.082516.tar.gz", exdir = "RawMatrix")2.2 表达矩阵构建
关键步骤代码示例：
library(data.table) library(dplyr) # 读取样本信息 sample_info <- fread("gdc_sample_sheet.2024-06-06.tsv") sample_info$Barcode <- substr(sample_info$`Sample ID`, 1, 15) # 样本筛选（01=肿瘤，11=正常） filtered_samples <- sample_info %>% filter(!duplicated(Barcode)) %>% filter(grepl("01$|11$", Barcode)) # 初始化表达矩阵 expr_matrix <- data.frame(gene_id=character(), gene_name=character(), gene_type=character())2.3 数据合并与过滤
合并所有样本的表达数据：
for (i in 1:nrow(filtered_samples)) { file_path <- paste0("RawMatrix/", filtered_samples$`File ID`[i], "/", filtered_samples$`File Name`[i]) sample_data <- fread(file_path) # 提取count数据 sample_counts <- sample_data[!1:4, c("gene_id", "unstranded")] colnames(sample_counts)[2] <- filtered_samples$Barcode[i] # 合并到主矩阵 if (i == 1) { expr_matrix <- sample_data[!1:4, c("gene_id", "gene_name", "gene_type")] } expr_matrix <- merge(expr_matrix, sample_counts, by="gene_id") } # 过滤低表达基因 expr_matrix <- expr_matrix[rowMeans(expr_matrix[, -c(1:3)]) > 1, ]3. 差异表达分析实战
3.1 使用edgeR进行差异分析
完整分析流程代码：
library(edgeR) # 准备分组信息 tumor_samples <- grep("-01$", colnames(expr_matrix), value=TRUE) normal_samples <- grep("-11$", colnames(expr_matrix), value=TRUE) group <- factor(c(rep("tumor", length(tumor_samples)), rep("normal", length(normal_samples)))) # 创建DGEList对象 dge <- DGEList(counts=expr_matrix[, c(tumor_samples, normal_samples)], genes=expr_matrix[, 1:3], group=group) # 过滤低表达基因 keep <- filterByExpr(dge) dge <- dge[keep, , keep.lib.sizes=FALSE] # 标准化 dge <- calcNormFactors(dge) # 差异分析 design <- model.matrix(~group) dge <- estimateDisp(dge, design) fit <- glmQLFit(dge, design) qlf <- glmQLFTest(fit) # 提取结果 deg_results <- topTags(qlf, n=Inf)$table3.2 结果解读与筛选
设置差异基因标准：
| 指标 | 阈值 | 生物学意义 |
|---|---|---|
| logFC | >2或<-2 | 表达量变化倍数 |
| FDR | <0.05 | 统计显著性 |
| 表达水平 | >100 | 确保有生物学意义 |
筛选差异基因代码：
deg_results$DEG <- "None" deg_results$DEG[deg_results$logFC > 2 & deg_results$FDR < 0.05] <- "Up" deg_results$DEG[deg_results$logFC < -2 & deg_results$FDR < 0.05] <- "Down" # 保存结果 write.csv(deg_results, "Results/LIHC_DEG_results.csv", row.names=FALSE)4. 高级可视化：火山图绘制
4.1 基础火山图实现
使用ggplot2绘制专业级火山图：
library(ggplot2) library(ggrepel) # 准备数据 deg_results$log10FDR <- -log10(deg_results$FDR) top_genes <- deg_results %>% group_by(DEG) %>% top_n(10, abs(logFC)) %>% pull(gene_name) # 绘制火山图 ggplot(deg_results, aes(x=logFC, y=log10FDR, color=DEG)) + geom_point(alpha=0.6, size=2) + scale_color_manual(values=c(Down="blue", None="gray", Up="red")) + geom_vline(xintercept=c(-2, 2), linetype="dashed") + geom_hline(yintercept=-log10(0.05), linetype="dashed") + geom_text_repel(data=subset(deg_results, gene_name %in% top_genes), aes(label=gene_name), size=3, box.padding=0.5) + labs(x="log2 Fold Change", y="-log10(FDR)", title="LIHC Tumor vs Normal Differential Expression") + theme_minimal() + theme(legend.position="bottom")4.2 图表美化技巧
提升出版级图表质量的几个关键参数：
# 高级定制化参数 volcano_plot <- volcano_plot + theme( plot.title = element_text(size=14, face="bold", hjust=0.5), axis.title = element_text(size=12), legend.text = element_text(size=10), panel.grid.major = element_line(color="gray90"), panel.background = element_rect(fill="white") ) + scale_x_continuous(breaks=seq(-10, 10, 2)) + coord_cartesian(xlim=c(-10, 10)) # 保存高质量图片 ggsave("Results/LIHC_volcano_plot.pdf", plot=volcano_plot, width=8, height=6, dpi=300)5. 常见问题排查与优化
5.1 报错解决方案
常见错误及解决方法：
包安装失败：
- 检查R版本是否过旧
- 尝试更换CRAN镜像源
- 对于Bioconductor包，确保使用BiocManager安装
内存不足：
# 增加内存限制 options(future.globals.maxSize=8000*1024^2)文件路径错误：
- 使用绝对路径替代相对路径
- 检查文件权限
- 确保文件名无特殊字符
5.2 性能优化建议
处理大型TCGA数据集时：
- 使用data.table替代data.frame提升读取速度
- 分批处理样本，避免内存溢出
- 考虑使用RDS格式保存中间结果
# 高效数据保存与读取 saveRDS(expr_matrix, "ProcessedData/expr_matrix.rds") expr_matrix <- readRDS("ProcessedData/expr_matrix.rds")5.3 分析流程自动化
将整个流程封装为函数：
run_deg_analysis <- function(project="LIHC", filter_expr=1, logfc_threshold=2, fdr_threshold=0.05) { # 包含所有分析步骤 # ... return(list(deg_results=deg_results, volcano_plot=volcano_plot)) } # 调用函数 lihc_results <- run_deg_analysis()