企业官网建设流程全解析

简介：这是一份聚焦ceRNA网络构建与分析的生物信息学报告，面向从事数据挖掘、肿瘤机制研究的初学者或科研人员，系统梳理了从TCGA数据下载、lncRNA/miRNA/mRNA差异表达到ceRNA调控网络构建的完整分析流程。资源整体为1个PDF文件，压缩包仅1.63MB，内容精炼，内含差异表达筛选表格（含logFC、logCPM、PValue、FDR）及网络分析关键步骤，便于快速掌握研究套路。目前已有85人学习下载。报告不仅覆盖热图绘制、靶向关系预测、ceRNA网络可视化，还进一步介绍生存曲线批量绘制与GO/KEGG功能富集分析，并结合差异基因筛选指标详解筛选方法。读者可据此理解癌症相关基因调控机制，也能将整套分析路径迁移至其他肿瘤数据中，适合作为生物信息入门及课题设计的参考模板。

1. 为什么 ceRNA 分析从 TCGA 差异表达开始

拿到一套 TCGA 的 RNA 测序数据，第一件事往往不是急着画网络图，而是先回答三个问题：lncRNA、miRNA、mRNA 这三类分子在肿瘤和正常样本里到底谁在显著变化。ceRNA（competing endogenous RNA）网络的核心假设是 lncRNA 通过竞争性结合 miRNA 来间接调控 mRNA 的表达，这个机制要成立，前提是参与网络的分子本身必须在表达层面有明显差异。换句话说，差异表达筛选是 ceRNA 分析的地基，地基没打牢，后面预测出来的调控关系大概率是噪声。

本文提到的分析报告覆盖了从 TCGA 数据下载、三种 RNA 差异表达、热图、靶向关系预测、ceRNA 网络构建、批量生存曲线到 GO/KEGG 富集分析的完整流程，非常适合正在做肿瘤生信入门、准备复现 ceRNA 文章套路的研究生和从业者。你会发现这套流程里最耗时间的不是绘图，而是筛选阈值的选择和数据格式的整理，下面按实际执行的顺序把每一步拆开讲。

2. TCGA 数据下载与差异表达三重筛选

2.1 用 TCGAbiolinks 下载三种 RNA 表达矩阵

TCGA 的数据分散在 GDC（Genomic Data Commons）门户里，lncRNA 和 mRNA 表达量来自 STAR-Counts 流程，miRNA 表达量来自 BCGSC miRNA 定量流程。手工在网页上勾选文件容易漏样本，而且下载下来的文件是 per-sample 格式，需要自己合并成矩阵，这一步用 R 的 TCGAbiolinks 包可以一次性完成。

library(TCGAbiolinks) # 下载 mRNA 和 lncRNA 的转录组表达矩阵 query_mrna <- GDCquery( project = "TCGA-COAD", data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts", sample.type = c("Primary Tumor", "Solid Tissue Normal") ) GDCdownload(query_mrna) mrna_expr <- GDCprepare(query_mrna) # 下载 miRNA 表达矩阵 query_mirna <- GDCquery( project = "TCGA-COAD", data.category = "Transcriptome Profiling", data.type = "miRNA Expression Quantification", workflow.type = "BCGSC miRNA Profiling", sample.type = c("Primary Tumor", "Solid Tissue Normal") ) GDCdownload(query_mirna) mirna_expr <- GDCprepare(query_mirna)

project参数按癌种填写，比如 COAD 是结肠癌，BRCA 是乳腺癌，LUAD 是肺腺癌。sample.type里的Primary Tumor和Solid Tissue Normal分别对应肿瘤和癌旁正常样本，这一项直接决定后续差异分析的比较组，填错了后面全盘皆输。下载完成后建议把mrna_expr和mirna_expr存成 RData 文件，因为 GDC 服务偶尔会断连，重新下载非常耗时。

TCGAbiolinks 返回的是 SummarizedExperiment 对象，表达矩阵可以用assay()提取，样本信息用colData()提取。mRNA 和 lncRNA 实际上在同一份转录组矩阵里，区别在于基因注释来源，Ensembl ID 以 ENSG 开头的是蛋白编码基因，以 ENSG 加_开头的长链非编码 RNA 在 Gencode 注释里有单独的分类标签，后续通过注释文件拆分即可。

2.2 差异表达分析的统计指标与筛选阈值

差异表达分析输出你表格里的四列：logFC、logCPM、PValue、FDR，这四个指标各有分工。logFC 是 log2 倍数变化，正值代表上调，负值代表下调；logCPM 是每百万reads中该基因的log2计数，衡量表达丰度，防止低表达基因干扰；PValue 来自统计检验；FDR 是多重检验校正后的 P 值，因为一次比较几万个基因，直接用原始 P 值会得到大量假阳性。

报告附件里 TCEAL6 的 logFC 是 -8.60，说明这个基因在肿瘤里表达量下降了约 2 的 8.6 次方，是极显著的下调基因。FDR 低至 8.87E-79，意味着这个差异基本不可能是抽样误差造成的。操作时建议先按 FDR 排序筛选，再根据研究目的决定是否收紧 logFC 阈值，如果后续要做 ceRNA 网络，阈值太严可能导致网络节点太少，网络图撑不起来。

2.3 edgeR 计算三张差异表达表

常用的差异表达工具有 edgeR、DESeq2、limma，其中 edgeR 对转录组数据稳定性好，计算速度快，适合处理 TCGA 这种大样本量数据。报告里出现 logCPM 这个指标，说明原分析很可能用的是 edgeR 或 limma-voom 的 CPM 体系。

library(edgeR) library(SummarizedExperiment) # 提取count矩阵，基因名作为行名 count_matrix <- assay(mrna_expr) gene_info <- rowData(mrna_expr) # 构建 DGEList 对象，group 标注肿瘤和正常 group <- ifelse(colData(mrna_expr)$definition == "Solid Tissue Normal", "normal", "tumor") y <- DGEList(counts = count_matrix, group = group) # 过滤低表达基因：在至少一组样本中 CPM > 1 keep <- filterByExpr(y, group = group) y <- y[keep, , keep.lib.sizes = FALSE] # 标准化：TMM 方法消除文库大小差异 y <- calcNormFactors(y) # 估计离散度并做准似然比检验 y <- estimateDisp(y) fit <- glmQLFit(y, design = model.matrix(~ group)) res <- glmQLFTest(fit, coef = 2) # 输出结果表：logFC、logCPM、PValue、FDR out <- topTags(res, n = Inf)$table write.csv(out, "mRNA_diff.csv")

filterByExpr在两组样本数接近时能自动决定过滤阈值，避免把只在少数样本里表达的基因纳入统计。calcNormFactors用的是 TMM 算法，它的作用是校正不同样本测序深度和RNA组成差异，这一步漏掉会导致 logFC 虚高。coef = 2指定比较第二个系数，也就是 tumor 和 normal 的差异，如果你的 group 因子顺序不一样，这里要改成对应的列号。

同样的代码可以分别跑 lncRNA 和 miRNA 矩阵，只要把输入换成拆出来的 lncRNA 子矩阵和 miRNA 表达矩阵。注意 miRNA 表达矩阵的列名和 mRNA 矩阵不一定能对上，合并前用intersect()和match()统一样本顺序。

注意：edgeR 的 CPM 和 logCPM 单位不同，CPM 是原始每百万计数，logCPM 是 log2(CPM + 2)，解读时不用纠结具体数值，只看相对高低即可。

2.4 热图绘制与样本分组核对

差异基因筛选出来后会绘制热图，报告里说热图可以放在正文也可以放附件，属于常规做法。热图的价值不在于好看，而在于验证差异基因是否真的能区分肿瘤和正常两组样本，如果聚类结果和样本分组混乱，说明数据质量有问题，需要回到上游检查。

library(pheatmap) # 取显著差异基因的表达量，做Z-score标准化 sig_genes <- rownames(out)[out$FDR < 0.05 & abs(out$logFC) > 2] mat <- assay(mrna_expr)[sig_genes, ] mat_z <- t(scale(t(mat))) # 绘制热图，按样本定义分组注释 annotation_col <- data.frame(group = group) rownames(annotation_col) <- colnames(mrna_expr) pheatmap(mat_z, scale = "none", annotation_col = annotation_col, show_rownames = FALSE, fontsize_row = 6, cluster_cols = TRUE, cluster_rows = TRUE, color = colorRampPalette(c("navy", "white", "firebrick3"))(100) )

t(scale(t(mat)))对每个基因跨样本做 Z-score 标准化，让上调基因显示为红色、下调基因为蓝色，否则表达量高的基因会盖掉低表达基因的颜色变化。热图的列聚类如果出现肿瘤样本和正常样本交叉混合，要立刻检查是否有样本标签错位或者批次效应。绘制 miRNA 热图时建议直接用 logCPM 值，因为 miRNA 表达量整体偏低，Z-score 容易放大噪声。

3. 靶向关系预测与 ceRNA 网络可视化

3.1 靶基因预测工具的组合策略

ceRNA 网络需要两类靶向关系：miRNA 靶向 mRNA 和 miRNA 靶向 lncRNA。miRNA 对 mRNA 的靶向关系预测工具很成熟，TargetScan 依据种子区（seed region）序列互补和保守性打分，miRDB 用机器学习模型预测，两者结果合并取交集可以显著降低假阳性。lncRNA 层面的靶向关系支持相对弱一些，因为 lncRNA 数据库注释不全，常用的是 starBase/ENCORI 的实验验证数据和 miRcode 的预测数据。

操作上先把差异 miRNA 列表分别输入 TargetScan 和 miRDB，下载各自的靶基因表，用 R 取交集后得到候选 mRNA 靶基因，再和差异 mRNA 列表取交集，这时候筛出来的才是既被预测靶向、又在表达上发生变化的基因。lncRNA 这边用 starBase 的预测结果直接和差异 lncRNA 取交集，因为实验支持的结合位点更有说服力。

3.2 构建 ceRNA 网络的核心逻辑：共享 miRNA 响应元件

ceRNA 机制的核心是 miRNA 响应元件（MRE），一个 mRNA 的 3'UTR 上有多个 miRNA 结合位点，lncRNA 如果也含有这些位点，就能竞争性吸附 miRNA，从而释放 mRNA。因此网络构建并不是把所有靶向关系都连上，而是要找到 lncRNA 和 mRNA 共享同一批 miRNA 的三角关系。例如 hsa-mir-145 同时靶向 mRNA A 和 lncRNA B，A 和 B 之间就构成 ceRNA 关系。

实际操作时用 R 做三元组匹配，以 miRNA 为桥梁把 lncRNA 和 mRNA 关联起来。差异 miRNA 列表里 hsa-mir-145、hsa-mir-133a、hsa-mir-1 这些肌肉相关 miRNA 在多个癌种中出现，意味着它们的靶基因和海绵 lncRNA 在肿瘤中往往协同变化。这一步生成的边文件是三列的 CSV，第一列是 lncRNA，第二列是 miRNA，第三列是 mRNA，Cytoscape 导入后会自动生成完整网络。

# 假设有三张表：mirna_mrna、mirna_lncrna、diff_mrna、diff_lncrna # 合并靶向关系并生成网络边表 library(dplyr) edges <- mirna_mrna %>% inner_join(mirna_lncrna, by = "miRNA") %>% filter(mRNA %in% diff_mrna$gene, lncRNA %in% diff_lncrna$gene) %>% select(lncRNA, miRNA, mRNA) write.csv(edges, "ceRNA_edges.csv", row.names = FALSE)

inner_join按 miRNA 列匹配两边的靶向关系，这一步会把没有同时出现在两张表里的 miRNA 自动剔除。过滤掉非差异基因后，网络规模会大幅缩小，留下的边都有表达差异和靶向预测的双重证据。生成的 CSV 建议再加一列interaction_type，值为lncRNA-miRNA或miRNA-mRNA，方便在 Cytoscape 里按关系类型设置线条颜色。

提示：如果发现网络边数太少，通常不是预测工具的问题，而是差异表达的阈值卡太严。把 FDR 从 0.05 放宽到 0.1，网络规模可能增加 30% 以上。

3.3 Cytoscape 网络可视化与核心节点识别

3.3.1 网络文件格式与导入

Cytoscape 导入边表时要求至少有 source node 和 target node 两列，后续的 miRNA 和 mRNA 类型可以用节点属性表单独维护。实际操作时我会准备两个文件，边的 CSV 里加入权重列，权重可以设为结合位点数量或预测得分，用于控制线条粗细。

# Cytoscape 导入步骤 # 1. File -> Import -> Network from File，选择 ceRNA_edges.csv # 2. 确认 Source Node 为 lncRNA 列，Target Node 为 miRNA 列 # 3. 将 miRNA 节点映射为三角形，lncRNA 映射为菱形，mRNA 映射为圆形 # 4. Style -> Fill Color 按节点类型分组着色 # 5. Layout -> yFiles Organic Layout，适合中等规模网络

网络里 mRNA 节点通常最多，如果图上 mRNA 把 miRNA 和 lncRNA 都盖住了，可以在样式面板里调小 mRNA 节点尺寸，并降低透明度。边样式里把 miRNA-mRNA 边设为虚线，lncRNA-miRNA 边设为实线，视觉上能快速区分两种调控关系。

3.3.2 从网络图中提取核心节点

网络建好后不能只截图发文章，还要用网络拓扑指标说明哪些节点是关键调控因子。Cytoscape 内置的 NetworkAnalyzer 可以计算度（degree）、介数中心性（betweenness centrality）、紧密度（closeness centrality），度高的 miRNA 说明它被大量 mRNA 和 lncRNA 竞争性结合，处于调控枢纽位置；介数中心性高的节点说明它连接着多个功能模块，删除它网络会分裂成多个小社区。

实际操作时我会把 NetworkAnalyzer 的结果导出，按 degree 排序选出前 10 个节点，再结合差异表达倍数筛选核心候选分子。报告中的 hsa-mir-10b、hsa-mir-145、hsa-mir-143 这类差异倍数高且已知有肿瘤抑制功能的 miRNA，如果在网络里度也排在前列，就可以作为后续生存分析和功能富集的重点对象。这里需要留意 miRNA 命名中的 hsa-mir 与 hsa-miR 区别，前者表示前体 miRNA，网络构建时应统一为成熟体命名 hsa-miR-145，否则同一个 miRNA 会被拆成多个节点，导致网络边数虚高。

4. 批量生存分析与 GO/KEGG 富集实操

4.1 批量生存曲线的 R 语言实现

ceRNA 网络里的节点必须和临床预后挂钩才有文章价值，这一步要做的是对网络里的每个 lncRNA、miRNA、mRNA 分组做生存曲线。TCGA 的临床信息里有 OS（总生存期）和 OS.time，通常取中位数把样本分成高表达组和低表达组，用 log-rank 检验比较两组生存差异。批量分析时用循环把每个基因跑一遍，把显著结果保存到汇总表。

library(survival) library(survminer) # 假设 expr_matrix 的行是基因，列是样本；clin 有 OS、OS.time genes <- rownames(expr_matrix) result <- data.frame() for (g in genes) { expr <- as.numeric(expr_matrix[g, ]) group <- ifelse(expr > median(expr), "high", "low") tmp <- data.frame(time = clin$OS.time, status = clin$OS, group = group) fit <- survfit(Surv(time, status) ~ group, data = tmp) p <- surv_pvalue(fit)$pval # 只保留显著结果的生存曲线，文件名以基因名命名 if (p < 0.05) { pdf(paste0("KM_", g, ".pdf"), width = 5, height = 5) print(ggsurvplot(fit, pval = TRUE, risk.table = TRUE)) dev.off() } result <- rbind(result, data.frame(gene = g, pvalue = p)) } result$FDR <- p.adjust(result$pvalue, method = "BH") write.csv(result, "survival_summary.csv")

median()是按基因表达量中位数分组的默认做法，但有些基因表达量偏态分布，中位数分组会让两组样本量严重失衡，建议用surv_cutpoint()函数自动寻找最优分组截点。surv_pvalue()计算的是 log-rank 检验 P 值，如果样本量小于 30，P 值会偏低，文章里一般会同时报告 Cox 回归的 HR 值。批量输出的 PDF 数量可能达到几百个，建议只在 P 值小于 0.05 时输出文件，其他情况只记录数值，避免生成大量无用的空图。

4.2 GO 富集与 KEGG 通路分析的参数设置

差异基因的功能注释用 clusterProfiler 完成，GO 富集覆盖生物过程（BP）、细胞组分（CC）、分子功能（MF）三个维度，KEGG 富集则关联代谢和信号通路。做 ceRNA 分析时建议分别对 mRNA 和 lncRNA 的靶基因做富集，因为 lncRNA 本身不编码蛋白，直接对 lncRNA 做 GO 没有意义，需要先通过 ceRNA 关系映射到它调控的 mRNA 靶基因。

library(clusterProfiler) library(org.Hs.eg.db) # 把差异mRNA的Symbol转成Entrez ID genes_entrez <- bitr(diff_mrna$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db ) # GO 富集 go_res <- enrichGO( gene = genes_entrez$ENTREZID, OrgDb = org.Hs.eg.db, ont = "BP", # BP/CC/MF 三选一 pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE ) # KEGG 富集：需要联网访问 KEGG API kegg_res <- enrichKEGG( gene = genes_entrez$ENTREZID, organism = "hsa", pvalueCutoff = 0.05, qvalueCutoff = 0.05 ) write.csv(as.data.frame(go_res), "GO_BP.csv") write.csv(as.data.frame(kegg_res), "KEGG.csv")

ont参数决定富集的 GO 子本体，BP 的注释条目最全，CC 和 MF 结果通常偏少，建议三个子本体都跑一遍，文章中优先展示 BP。readable = TRUE把 Entrez ID 转回基因 Symbol，否则表格里全是数字 ID，没法直接阅读。KEGG 富集依赖 KEGG 服务器，国内网络环境下经常连接超时，可以用download.kegg()将 pathway 数据缓存到本地，或者改用 ReactomePA 包做通路富集作为替代。

4.3 富集结果的可视化与结果解读

富集分析跑完以后，需要把 GO 和 KEGG 的结果画成点图和网络图。点图是文章的标配，横轴是 GeneRatio，纵轴是通路名称，点大小代表富集到的基因数，颜色代表 P 值。网络图能看出哪些基因同时富集到多条通路，也就是多效基因，这类基因往往是 ceRNA 网络里值得重点关注的下游效应分子。

library(ggplot2) # GO 点图 dotplot(go_res, showCategory = 15, font.size = 10) + ggtitle("GO Biological Process Enrichment") # KEGG 通路关系网络图 cnetplot(kegg_res, showCategory = 5, foldChange = diff_mrna$logFC)

showCategory控制图中显示的通路数量，建议保留 10 到 15 条，太多图会变成一团黑。cnetplot的foldChange参数传入 logFC 向量后，节点颜色会按表达变化方向渐变，图中就能看出哪些通路的靶基因整体上调、哪些整体下调。由于富集通路较多，文章里通常只保留和癌症相关的通路，比如 p53 信号通路、细胞凋亡、PI3K-Akt 信号通路等，不用把所有显著通路都列出来。报告正文里强调的 GO 和 KEGG 是 ceRNA 分析的收尾步骤，它的作用是验证网络中的分子是否富集在已知的癌症通路里，如果富集结果全是核糖体这类管家通路，就要考虑差异基因筛选是否过宽。

5. ceRNA 分析排错与批量提速的实用技巧

5.1 数据批次效应与分组不一致的排查

TCGA 数据下载后最常踩的坑是样本类型编码问题，colData(mrna_expr)$definition里除了Primary Tumor和Solid Tissue Normal，还可能出现Metastatic或者Additional New Primary，这些样本如果不剔除会被错误归入肿瘤组，直接污染差异表达结果。执行差异分析前先检查分组计数。

table(colData(mrna_expr)$definition)

如果结果里出现第三类及以上样本类型，用colData(mrna_expr)$definition %in% c("Primary Tumor", "Solid Tissue Normal")过滤后再构建 DGEList。另外，同一批 TCGA 项目里不同平台测序的数据不能直接合并，GDC 查询时要注意platform参数是否一致，Illumina HiSeq 和 NovaSeq 的数据混用会引入明显的批次效应，表现为热图聚类先按平台分开而不是按肿瘤正常分开。

5.2 差异基因太少时的调整策略

严格按照 FDR < 0.05 且 |logFC| > 2 筛选，某些癌种可能只得到几十个差异基因，这样的基因量构建 ceRNA 网络后只有 20 个节点，没有分析价值。我一般会分两步调整：先检查 logCPM 过滤是否过严，filterByExpr在样本量差异大时会把只在正常组高表达的基因过滤掉，可以把过滤条件改为rowSums(cpm(y) > 1) >= 2，保留更多低丰度基因；如果还是太少，把 |logFC| 阈值放宽到 1，FDR 放宽到 0.1，同时在文章中注明筛选标准是 FDR < 0.1。

还有一种情况是差异 miRNA 数量少。TCGA 的 miRNA 定量结果是 counts 矩阵，但很多 miRNA 在肿瘤和正常组织里表达都非常低，导致检验功效不足。这时可以改用 miRNA 的 RPM（reads per million）值做 limma 差异分析，limma 对低表达基因的 P 值更为宽松，实测下来 miRNA 差异基因数量通常能提升不少。

5.3 批量生存分析与数据缓存优化

批量生存分析跑几百个基因时，主要瓶颈不在统计检验，而在反复读取表达矩阵和磁盘 IO。我的做法是用 data.table 的fread把表达矩阵一次性读入内存，转成矩阵后每次按行索引取值，避免循环里多次读 CSV 文件。代码里改成expr_matrix[g, ]直接索引比expr_matrix[g, , drop = FALSE]速度更快，因为后者会创建不必要的副本。

library(data.table) # 用 fread 加速读取，并转成矩阵 expr_matrix <- as.matrix(fread("expr_matrix.csv", data.table = FALSE)) rownames(expr_matrix) <- expr_matrix[, 1] expr_matrix <- expr_matrix[, -1]

对于特别大的矩阵，建议把差异基因列表先写好，只对差异基因做批量生存分析，不用对全部两万个基因都跑一遍。p 值校正时用p.adjust(pvalue, method = "BH")，校正后仍然显著的基因数量会明显下降，文章中只需报告这些关键基因的生存曲线即可。最后的生存汇总表记得加上基因类型列，比如 lncRNA、miRNA、mRNA，方便后期按类别统计显著比例。

本文还有配套的精品资源，点击获取