R语言实战：从TCGA官网下载到火山图，手把手搞定肝癌(LIHC)差异表达分析全流程

R语言实战从TCGA官网下载到火山图手把手搞定肝癌(LIHC)差异表达分析全流程在生物信息学研究中TCGA数据库是癌症基因组分析的重要资源宝库。对于肝癌(LIHC)研究而言掌握从原始数据获取到差异表达分析的全流程是每个科研人员必备的核心技能。本文将带你用R语言完成从TCGA数据下载、预处理、差异分析到可视化呈现的完整过程特别针对初学者可能遇到的各类问题提供解决方案。1. TCGA数据获取与前期准备1.1 数据下载流程详解访问TCGA官方数据门户(https://portal.gdc.cancer.gov/)后按以下步骤操作在Cohort Builder中选择Program为TCGAProject为LIHC进入Repository后在侧边栏设置筛选条件Experimental Strategy: RNA-SeqData Category: Transcriptome profilingData 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_mirrorhttps://mirrors.tuna.tsinghua.edu.cn/bioconductor)包依赖冲突时建议新建干净的R session2. 数据整理与质量控制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_idcharacter(), gene_namecharacter(), gene_typecharacter())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, bygene_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), valueTRUE) normal_samples - grep(-11$, colnames(expr_matrix), valueTRUE) group - factor(c(rep(tumor, length(tumor_samples)), rep(normal, length(normal_samples)))) # 创建DGEList对象 dge - DGEList(countsexpr_matrix[, c(tumor_samples, normal_samples)], genesexpr_matrix[, 1:3], groupgroup) # 过滤低表达基因 keep - filterByExpr(dge) dge - dge[keep, , keep.lib.sizesFALSE] # 标准化 dge - calcNormFactors(dge) # 差异分析 design - model.matrix(~group) dge - estimateDisp(dge, design) fit - glmQLFit(dge, design) qlf - glmQLFTest(fit) # 提取结果 deg_results - topTags(qlf, nInf)$table3.2 结果解读与筛选设置差异基因标准指标阈值生物学意义logFC2或-2表达量变化倍数FDR0.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.namesFALSE)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(xlogFC, ylog10FDR, colorDEG)) geom_point(alpha0.6, size2) scale_color_manual(valuesc(Downblue, Nonegray, Upred)) geom_vline(xinterceptc(-2, 2), linetypedashed) geom_hline(yintercept-log10(0.05), linetypedashed) geom_text_repel(datasubset(deg_results, gene_name %in% top_genes), aes(labelgene_name), size3, box.padding0.5) labs(xlog2 Fold Change, y-log10(FDR), titleLIHC Tumor vs Normal Differential Expression) theme_minimal() theme(legend.positionbottom)4.2 图表美化技巧提升出版级图表质量的几个关键参数# 高级定制化参数 volcano_plot - volcano_plot theme( plot.title element_text(size14, facebold, hjust0.5), axis.title element_text(size12), legend.text element_text(size10), panel.grid.major element_line(colorgray90), panel.background element_rect(fillwhite) ) scale_x_continuous(breaksseq(-10, 10, 2)) coord_cartesian(xlimc(-10, 10)) # 保存高质量图片 ggsave(Results/LIHC_volcano_plot.pdf, plotvolcano_plot, width8, height6, dpi300)5. 常见问题排查与优化5.1 报错解决方案常见错误及解决方法包安装失败检查R版本是否过旧尝试更换CRAN镜像源对于Bioconductor包确保使用BiocManager安装内存不足# 增加内存限制 options(future.globals.maxSize8000*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(projectLIHC, filter_expr1, logfc_threshold2, fdr_threshold0.05) { # 包含所有分析步骤 # ... return(list(deg_resultsdeg_results, volcano_plotvolcano_plot)) } # 调用函数 lihc_results - run_deg_analysis()