基于 R 语言的 miRNA 数据分析：靶基因预测与通路富集探究

基于 R 语言的 miRNA 数据分析：靶基因预测与通路富集探究
一、引言
在生命科学领域，微小核糖核酸（miRNA）作为一类长度较短的非编码 RNA 分子，近年来备受关注。miRNA 能够通过与靶信使核糖核酸（mRNA）的互补配对，在转录后水平对基因表达进行精细调控，进而影响细胞的增殖、分化、凋亡等多种生理过程，与众多疾病的发生发展密切相关 。准确识别 miRNA 的靶基因以及深入探究其参与的生物学通路，对于揭示 miRNA 的功能机制至关重要。随着生物信息学技术的飞速发展，利用 R 语言这一强大的数据分析工具进行 miRNA 相关研究成为主流趋势。R 语言拥有丰富的软件包，能够高效地处理和分析海量的生物数据，为 miRNA 研究提供了有力支持。本研究旨在运用 R 语言，借助 multiMiR 等工具包，对特定 miRNA 的靶基因进行预测，并通过 KEGG 和 GO 通路富集分析，深入探究这些 miRNA 在生物学过程中的潜在作用机制，为相关领域的研究提供有价值的参考。
二、材料与方法
2.1 实验材料
本研究选取了 19 种人类 miRNA，具体包括 “hsa - miR - 9 - 5p”、“hsa - miR - 21 - 5p”、“hsa - miR - 215 - 5p” 等。这些 miRNA 是基于前期研究报道以及在相关生物学过程中可能发挥关键作用而选定。
2.2 实验方法
2.2.1 数据处理与分析工具
使用 R 语言（版本 X.X.X）进行数据处理和分析，主要借助以下工具包：
r
library(multiMiR)
library(clusterProfiler)
library(org.Hs.eg.db)
library(dplyr)
在开始分析前，确保所有相关包均已成功安装并加载至 R 环境中。
2.2.2 miRNA 靶基因预测
将选定的 19 种 miRNA 构建成列表对象mir_list_converted：
r
mir_list_converted <- list(
“hsa - miR - 9 - 5p”,
“hsa - miR - 21 - 5p”,
“hsa - miR - 215 - 5p”,
“hsa - miR - 222 - 5p”,
“hsa - miR - 221 - 5p”,
“hsa - miR - 10b - 5p”,
“hsa - miR - 26a - 5p”,
“hsa - miR - 340 - 5p”,
“hsa - miR - 551 - 5p”,
“hsa - miR - 17 - 5p”,
“hsa - miR - 210 - 5p”,
“hsa - miR - 34a - 5p”,
“hsa - miR - 125b - 5p”,
“hsa - miR - 141 - 5p”,
“hsa - miR - 124 - 5p”,
“hsa - miR - 203 - 5p”,
“hsa - miR - 200 - 5p”,
“hsa - miR - 223 - 5p”,
“hsa - miR - 132 - 5p”
)
通过循环遍历该列表，利用get_multimir函数从 miRDB 数据库中查询每个 miRNA 的靶基因。该函数设置参数org = 'hsa’指定物种为人类，table = 'all’表示获取所有相关数据：
r

创建一个空列表来存储结果

mir_targets <- list()

遍历miRNA列表

for (mir in mir_list_converted) {

查询miRNA的靶基因

result <- get_multimir(mirna = mir, org = ‘hsa’, table = ‘all’)

将结果存储到列表中

mir_targets[[mir]] <- result@data
}
为初步了解数据情况，使用head函数查看了mir_targets[[“miR - 9 - 5p”]]的前几行数据，通过sapply函数统计了每个 miRNA 靶基因的数量 ：
r

查看某个miRNA的靶基因

head(mir_targets[[“hsa - miR - 9 - 5p”]])

查看所有miRNA的靶基因数量

sapply(mir_targets, nrow)
进一步将所有 miRNA 靶基因的符号提取出来，构建成数据框df，并导出为 CSV 文件all_target_symbols.csv，以便后续分析。同时，计算所有 miRNA 靶基因的交集，得到交集基因列表intersected_genes，这有助于筛选出可能受多种 miRNA 共同调控的关键基因：
r
all_target_symbols <- lapply(mir_targets, function(x) x$target_symbol)

将 all_target_symbols 转换为数据框

df <- data.frame(Target_Symbols = unlist(all_target_symbols))

导出为 CSV 文件

write.csv(df, file = “all_target_symbols.csv”, row.names = FALSE)

计算所有 miRNA 靶基因的交集

intersected_genes <- Reduce(intersect, all_target_symbols)
2.2.3 通路富集分析
首先，将所有 miRNA 靶基因符号合并成一个向量combined_target_symbols，去除其中的重复符号，得到唯一的基因符号列表unique_target_symbols。然后，从org.Hs.eg.db包中获取有效的人类基因符号valid_symbols，筛选出在unique_target_symbols中存在的有效目标基因符号valid_target_symbols。利用bitr函数将这些有效的基因符号转换为 Entrez ID，为后续的通路富集分析做准备 ：
r
combined_target_symbols <- unlist(all_target_symbols)

去除重复的基因符号

unique_target_symbols <- unique(combined_target_symbols)

获取 org.Hs.eg.db 中有效的基因符号

valid_symbols <- keys(org.Hs.eg.db, keytype = “SYMBOL”)

筛选出有效的目标基因符号

valid_target_symbols <- unique_target_symbols[unique_target_symbols %in% valid_symbols]

将有效的基因符号转换为 Entrez ID

entrez_ids <- bitr(valid_target_symbols, fromType = “SYMBOL”, toType = “ENTREZID”, OrgDb = org.Hs.eg.db)$ENTREZID
定义感兴趣的信号通路，包括 “TGF - beta.*signaling pathway”、“JAK - STAT.*signaling pathway”、“Notch.*signaling pathway”、“HIF - 1.*signaling pathway”、“NF - kappa B.*signaling pathway” 等：
r
pathways <- c(“TGF - beta.*signaling pathway”, “JAK - STAT.*signaling pathway”, “Notch.*signaling pathway”,
“HIF - 1.*signaling pathway”, “NF - kappa B.*signaling pathway”)
使用enrichKEGG函数对转换后的 Entrez ID 进行 KEGG 通路富集分析，设置参数organism = "hsa"指定物种为人类，pAdjustMethod = "BH"采用 Benjamini - Hochberg 方法进行多重检验校正，qvalueCutoff = 0.05设定校正后 p 值的阈值为 0.05 。将 KEGG 通路富集分析的结果导出为 CSV 文件kegg_enrichment_results.csv：
r

进行 KEGG 通路富集分析

kegg_enrich <- enrichKEGG(gene = entrez_ids, organism = “hsa”, pAdjustMethod = “BH”, qvalueCutoff = 0.05)
write.csv(kegg_enrich@result, file = “kegg_enrichment_results.csv”, row.names = FALSE)
同样地，利用enrichGO函数进行 GO 通路富集分析，设置参数OrgDb = org.Hs.eg.db指定使用的物种注释数据库，ont = "ALL"表示同时分析生物过程（BP）、分子功能（MF）和细胞组分（CC）三个方面，pAdjustMethod = "BH"和qvalueCutoff = 0.05与 KEGG 分析中的设置一致 。将 GO 通路富集分析的结果保存为 CSV 文件go_enrichment_results.csv：
r

进行 GO 通路富集分析

go_enrich <- enrichGO(gene = entrez_ids,
OrgDb = org.Hs.eg.db,
ont = “ALL”, # 可以选择 “BP”（生物过程）、“MF”（分子功能）、“CC”（细胞组分）或 “ALL”
pAdjustMethod = “BH”,
qvalueCutoff = 0.05)

保存 GO 富集分析结果到 CSV 文件

write.csv(go_enrich@result, file = “go_enrichment_results.csv”, row.names = FALSE)
为了聚焦于感兴趣的通路，使用grep函数从 KEGG 富集结果中筛选出与定义的通路名称匹配的条目，构建selected_pathways数据集。为直观展示富集分析结果，使用barplot函数绘制 KEGG 富集结果的条形图，展示前 20 个富集通路，以校正后的 p 值进行上色；使用dotplot函数绘制selected_pathways的气泡图，展示前 10 个富集通路，同样按校正后的 p 值上色，气泡大小表示基因比例 ：
r

筛选出感兴趣的通路结果

selected_pathways <- kegg_enrich[grep(paste(pathways, collapse = “|”), kegg_enrich$Description, ignore.case = TRUE), ]

绘制条形图

barplot(kegg_enrich,
showCategory = 20, # 显示前 20 个富集通路
title = “KEGG Pathway Enrichment Analysis”,
color = “p.adjust”) # 按校正后的 p 值上色

绘制气泡图

dotplot(selected_pathways,
showCategory = 10, # 显示前 10 个富集通路
title = “KEGG Pathway Enrichment Analysis”,
color = “p.adjust”, # 按校正后的 p 值上色
size = “geneRatio”) # 气泡大小表示基因比例