基因共表达网络分析的R包BioNERO

参考链接
https://bioconductor.org/packages//release/bioc/vignettes/BioNERO/inst/doc/vignette_01_GCN_inference.html

我理解的是 底层代码还是WGCNA的代码，但是进行了封装，使用起来更简单。增加了数据过滤的函数，过滤数据更容易。画图代码进行了更行，不再是WGCNA包的配色，出的结果更美观。

分析代码
教程里用到的R版本是4.5 BioNERO版本是1.18，我用的R版本是 4.2.3，用biocManager安装默认的版本是 1.6.1。这个版本按照上面的教程的代码稍微不太一样，出的图的配色还是WGCNA的经典配色

输入数据是一个SummarizedExperiment对象，我们自己转录组分析通常拿到的是一个基因表达矩阵，我们可以把示例数据里的基因表达矩阵提取出来，自己试着构建一个SummarizedExperiment对象

zma.se@assays@data@listData[[1]] %>%
rownames_to_column(“geneid”) %>%
write_tsv(“wgcna_example_exp.tsv”)

zma.se@colData %>% as.data.frame() %>%
rownames_to_column(“sampleid”) %>%
write_tsv(“wgcna_example_metadata.tsv”)

表达的数据格式是 行是基因，列是样本，第一列是geneid

meta 数据 每行是一个样本，列表示分组，这个示例数据分组是每个样本来源的组织 第一列是样本id

meta.dat<-read_tsv(“wgcna_example_metadata.tsv”) %>%
column_to_rownames(“sampleid”)

exp.dat<-read_tsv(“wgcna_example_exp.tsv”) %>%
column_to_rownames(“geneid”) %>%
dplyr::select(all_of(rownames(meta.dat)))

identical(rownames(meta.dat),colnames(exp.dat))

library(SummarizedExperiment)
se.example <- SummarizedExperiment(
assays = list(counts = exp.dat),
colData = meta.dat
)
首先看下数据里有没有缺失值

sum(is.na(se.example))
有缺失值的话可以给数据做个填充，这里默认的是填充为0

exp_filt <- replace_na(se.example)
表达矩阵数据通常很少会出现缺失值的情况

根据表达量对数据进行过滤

exp_filt <- remove_nonexp(se.example, method = “median”, min_exp = 10)
这里采取的过滤标准是 基因表达量中位数大于等于10

还可以根据基因表达的变异程度对数据进行过滤

exp_filt <- filter_by_variance(exp_filt, n = 2000)
过滤离群样本

exp_filt <- ZKfiltering(exp_filt, cor_method = “pearson”)
校正/调整混杂因素

exp_filt <- PC_correction(exp_filt)
我这里为了加快运行速度根据基因表达变异程度过滤，只保留变异排名前2000的基因

exp_filt <- filter_by_variance(se.example, n = 2000)
exp_filt <- ZKfiltering(exp_filt, cor_method = “pearson”)
exp_filt <- PC_correction(exp_filt)
一些数据探索

基于样本聚类热图

基于基因表达聚类热图

PCA结果图

p <- plot_heatmap(final_exp, type = “samplecor”, show_rownames = FALSE)
p
p <- plot_heatmap(
final_exp[1:50, ], type = “expr”, show_rownames = FALSE, show_colnames = FALSE
)
p
plot_PCA(final_exp)
计算软阈值

sft <- SFT_fit(final_exp, net_type = “signed hybrid”, cor_method = “pearson”)
这里有一个 net_type参数，之前一直没搞懂是啥意思，问了问deepseek

构建共表达网络

power <- sft$power

net <- exp2gcn(
final_exp, net_type = “signed hybrid”, SFTpower = power,
cor_method = “pearson”
)

names(net)

几个结果图

plot_dendro_and_colors(net)
plot_eigengene_network(net)

plot_ngenes_per_module(net)
模块特征相关性

MEtrait <- module_trait_cor(exp = final_exp, MEs = net$MEs)
plot_module_trait_cor(MEtrait)
还可以做富集分析，这里需要用到interproscan的注释结果，注释结果的格式，两列

data(zma.interpro)

zma.interpro %>% dim()

zma.interpro %>% head()

之前看到有的论文里会做结构域的富集，那用这个函数应该可以做

interpro_enrichment <- module_enrichment(
net = net,
background_genes = rownames(final_exp),
annotation = zma.interpro
)
还有一个针对特定基因集做富集分析的函数

enrichment_analysis()

每个模块的hub基因

hubs <- get_hubs_gcn(final_exp, net)

每个基因属于哪个模块

net$genes_and_modules
某个特定模块的网络形式可视化

edges_filtered <- get_edge_list(
net, module = “red”,
filter = TRUE, method = “min_cor”, rcutoff = 0.7
)

plot_gcn(
edgelist_gcn = edges_filtered,
net = net,
color_by = “module”,
hubs = hubs
)

之前看到一个R包 ggnetview 可视化wgcna的结果很好看，抽时间看看这个的输出结果怎么和ggnetview的结果衔接

欢迎大家关注我的公众号

小明的数据分析笔记本