BioNERO：快速执行基因共表达网络推理分析

原创生信碱移生信碱移 2023-06-25 22:16

老铁快点击蓝字 关注起来

生信碱移

BioNERO

本文介绍了一个基因共表达网络推理的r包：BioNERO

GitHub项目地址：

https://github.com/almeidasilvaf/BioNERO

参考文档：

http://bioconductor.org/packages/release/bioc/vignettes/BioNERO/inst/doc/vignette_01_GCN_inference.html

http://bioconductor.org/packages/release/bioc/vignettes/BioNERO/inst/doc/vignette_02_GRN_inference.html

http://bioconductor.org/packages/release/bioc/vignettes/BioNERO/inst/doc/vignette_03_network_comparison.html

迄今为止，已经开发了几个软件包来从表达数据推断基因共表达网络（GCN）和基因调控网络（GRN），例如 WGCNA、CEMiTool、petal和minet。然而，它们都无法处理网络分析工作流程的所有方面，并且用户需要使用其他包来构建标准的分析管道。此外，网络推断需要扎实的线性代数和统计学背景，导致缺乏经验的研究人员很难正确预处理其表达数据并从推断的网络中提取具有生物学意义的信息。

在这里，Almeida-Silva等人展示了 「BioNERO」，这是一个 R/Bioconductor 软件包，它将基因网络推理工作流程的所有步骤集成在一个软件包中。BioNERO使用最先进的方法来预处理表达数据，帮助研究者从表达数据推断 GCN 和 GRN并分析基因网络以进行生物学解释，甚至比较物种内和物种间的网络。此外，BioNERO可用于探索蛋白质-蛋白质相互作用网络的拓扑特性，例如中心识别和群落检测。

以下内容将介绍使用该包进行WGCNA分析：

安装方式

BioNERO可以通过以下命令安装：

if(!requireNamespace('BiocManager', quietly = TRUE))
  install.packages('BiocManager')

BiocManager::install("BioNERO")

数据加载和预处理

在本示例中，将使用在 TPM 中标准化的玉米基因表达数据：zma.se，数据集存储为 SummarizedExperiment 对象。需要注意的是BioNERO输入的表达数据可以是 SummarizedExperiment 对象，也可以是基因表达矩阵或数据框，行为基因，列为样本。

data(zma.se)

# Take a quick look at the data
zma.se
## class: SummarizedExperiment 
## dim: 10802 28 
## metadata(0):
## assays(1): ''
## rownames(10802): ZeamMp030 ZeamMp044 ... Zm00001d054106 Zm00001d054107
## rowData names(0):
## colnames(28): SRX339756 SRX339757 ... SRX2792103 SRX2792104
## colData names(1): Tissue
SummarizedExperiment::colData(zma.se)
## DataFrame with 28 rows and 1 column
##                    Tissue
##               <character>
## SRX339756       endosperm
## SRX339757       endosperm
## SRX339758       endosperm
## SRX339762       endosperm
## SRX339763       endosperm
## ...                   ...
## SRX2792107 whole_seedling
## SRX2792108 whole_seedling
## SRX2792102 whole_seedling
## SRX2792103 whole_seedling
## SRX2792104 whole_seedling

自动、一步数据预处理

用户可以选择快速自动的数据预处理命令：

final_exp <- exp_preprocess(
    zma.se, min_exp = 10, variance_filter = TRUE, n = 2000
)
## Number of removed samples: 1
identical(dim(exp_filt), dim(final_exp))
## [1] TRUE

# Take a look at the final data
final_exp
## class: SummarizedExperiment 
## dim: 2000 27 
## metadata(0):
## assays(1): ''
## rownames(2000): ZeamMp030 ZeamMp092 ... Zm00001d054093 Zm00001d054107
## rowData names(0):
## colnames(27): SRX339756 SRX339757 ... SRX2792103 SRX2792104
## colData names(1): Tissue

注意：在这里，我们使用 TPM 标准化数据。如果您的表达式数据作为原始读取计数count，请在 exp_preprocess() 中设置参数 「vstransform = TRUE」。这会将 DESeq2 的方差稳定变换应用于您的count数据。

此外，用户可以也使用单个函数逐步预处理数据。上方函数exp_preprocess() 是函数replace_na()、remove_nonexp()、filter_by_variance()、ZKfiltering() 和PC_ Correction() 的集成。

更多详细请参考作者文档：

http://bioconductor.org/packages/release/bioc/vignettes/BioNERO/inst/doc/vignette_01_GCN_inference.html。

探索性数据分析

BioNERO包括一些方便数据探索的功能。创建这些函数是为了避免键入虽然很小但会多次使用的代码块。这里的想法是让生物网络分析的用户体验尽可能简单。

绘制热图：

函数plot_heatmap()在绘制样本或基因表达之间相关性的热图。

# Heatmap of sample correlations
p <- plot_heatmap(final_exp, type = "samplecor", show_rownames = FALSE)
p


# Heatmap of gene expression (here, only the first 50 genes)
p <- plot_heatmap(
    final_exp[1:50, ], type = "expr", show_rownames = FALSE, show_colnames = FALSE
)
p

主成分（PCA）分析

函数plot_PCA()执行 PCA 并绘制用户选择的任意一对 PC（默认情况下为 PC1 和 PC2）以及每个 PC 的解释方差百分比

plot_PCA(final_exp)

推断基因共表达网络

现在我们有了过滤和归一化的表达数据，我们可以使用 WGCNA 算法重建「基因共表达网络 (GCN」)。首先，我们需要确定最合适的β幂，使网络满足无尺度分布。我们使用函数 SFT_fit() 来做到这一点。相关性被提升到 β 次方以放大它们的距离，从而使模块检测算法更加准确。需要注意的是，β值越大，网络越接近无尺度分布。然而，非常高的 β 幂会降低平均连通性，这是我们不想看到的。为了解决这个权衡问题，先前的研究总是选择高于特定阈值（默认为 0.8）的最低 β 幂。这使得网络接近无尺度分布，而不会显著降低平均连通性。

sft <- SFT_fit(final_exp, net_type = "signed hybrid", cor_method = "pearson")
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      3    0.220 -0.218          0.178   278.0    303.00  598.0
## 2      4    0.416 -0.382          0.272   196.0    199.00  472.0
## 3      5    0.573 -0.468          0.462   145.0    136.00  381.0
## 4      6    0.675 -0.536          0.584   110.0     95.70  312.0
## 5      7    0.748 -0.584          0.676    86.3     70.00  259.0
## 6      8    0.791 -0.653          0.735    68.8     51.90  221.0
## 7      9    0.803 -0.717          0.761    55.8     38.60  191.0
## 8     10    0.815 -0.775          0.790    45.8     29.90  167.0
## 9     11    0.821 -0.828          0.815    38.1     22.90  147.0
## 10    12    0.838 -0.874          0.850    32.0     17.90  130.0
## 11    13    0.847 -0.913          0.876    27.2     14.30  116.0
## 12    14    0.856 -0.943          0.893    23.2     11.80  104.0
## 13    15    0.875 -0.973          0.913    20.0      9.79   93.0
## 14    16    0.892 -0.997          0.937    17.3      8.00   83.9
## 15    17    0.897 -1.020          0.941    15.1      6.75   76.0
## 16    18    0.891 -1.070          0.948    13.3      5.79   69.7
## 17    19    0.888 -1.100          0.950    11.7      4.96   64.2
## 18    20    0.888 -1.130          0.957    10.4      4.27   59.4
sft$power
## [1] 9
power <- sft$power
sft$plot

现在，我们可以使用计算出的β幂来推断 GCN。运行函数exp2gcn()推断 GCN 并输出 7 个元素的列表，每个元素都将被分析管道中的其他函数使用。

net <- exp2gcn(
    final_exp, net_type = "signed hybrid", SFTpower = power, 
    cor_method = "pearson"
)
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
names(net)
## [1] "adjacency_matrix"    "MEs"                 "genes_and_modules"  
## [4] "kIN"                 "correlation_matrix"  "params"             
## [7] "dendro_plot_objects"

函数exp2gcn()将对象保存在结果列表的最后一个元素中，随后可用于绘制 GCN 论文中的常见图形。这些图形符合出版要求。

# Dendro and colors
plot_dendro_and_colors(net)

# Eigengene networks
plot_eigengene_network(net)

让我们看看每个模块的基因数量

plot_ngenes_per_module(net)

探索基因共表达网络

评估模块稳定性

module_stability(final_exp, net, nRuns = 5)
##  ...working on run 1 ..
##  ...working on run 2 ..
##  ...working on run 3 ..
##  ...working on run 4 ..
##  ...working on run 5 ..
##  ...working on run 6 ..

模块-性状关联分析

MEtrait <- module_trait_cor(exp = final_exp, MEs = net$MEs)
head(MEtrait)
##        ME          trait          cor    pvalue  group
## 1 MEblack      endosperm -0.166480994 0.4065715 Tissue
## 2 MEblack         pollen  0.213691004 0.2845053 Tissue
## 3 MEblack whole_seedling -0.020227505 0.9202318 Tissue
## 4 MEbrown      endosperm  0.003843583 0.9848197 Tissue
## 5 MEbrown         pollen -0.020729547 0.9182584 Tissue
## 6 MEbrown whole_seedling  0.012815311 0.9494147 Tissue
plot_module_trait_cor(MEtrait)

可视化模块表达谱

plot_expression_profile(
    exp = final_exp, 
    net = net, 
    plot_module = TRUE, 
    modulename = "yellow"
)

富集分析

# Enrichment analysis for conserved protein domains (Interpro)
data(zma.interpro)
interpro_enrichment <- module_enrichment(
    net = net, 
    background_genes = rownames(final_exp),
    annotation = zma.interpro
)
## Enrichment analysis for module black...
## Enrichment analysis for module brown...
## Enrichment analysis for module darkgreen...
## Enrichment analysis for module darkolivegreen...
## Enrichment analysis for module greenyellow...
## Enrichment analysis for module lightyellow...
## Enrichment analysis for module midnightblue...
## Enrichment analysis for module paleturquoise...
## Enrichment analysis for module violet...
## Enrichment analysis for module yellow...

# Print results without geneIDs for better visualization
interpro_enrichment[, -6]
##                                        term genes all         pval         padj
## 185                      Histone H2A/H2B/H3    43  44 2.155952e-09 4.840112e-07
## 186                             Histone H2B    14  14 6.217795e-04 3.489738e-02
## 187                       Histone H3/CENP-A    15  15 3.659921e-04 2.347578e-02
## 188                              Histone H4    15  15 3.659921e-04 2.347578e-02
## 189                            Histone-fold    58  60 1.083394e-11 4.864438e-09
## 330           Ribosomal protein L2 domain 2    18  18 7.448332e-05 6.688602e-03
## 395     Translation protein SH3-like domain    22  22 8.872064e-06 1.327852e-03
## 396 Translation protein, beta-barrel domain    26  27 1.235834e-05 1.387224e-03
## 301                   Protein kinase domain     5  18 1.202246e-04 2.699043e-02
## 446         Zinc finger, RING/FYVE/PHD-type     5  18 1.202246e-04 2.699043e-02
## 53                    Aquaporin transporter     3   5 9.644015e-05 4.330163e-02
##           module
## 185        black
## 186        black
## 187        black
## 188        black
## 189        black
## 330        black
## 395        black
## 396        black
## 301  lightyellow
## 446  lightyellow
## 53  midnightblue

枢纽基因鉴定

hubs <- get_hubs_gcn(final_exp, net)
head(hubs)
##             Gene Module  kWithin
## 1 Zm00001d033147  black 188.3864
## 2 Zm00001d049790  black 181.4522
## 3 Zm00001d005649  black 180.7062
## 4 Zm00001d045448  black 180.6744
## 5 Zm00001d008203  black 178.7147
## 6 Zm00001d023340  black 177.7553

提取子模块

edges <- get_edge_list(net, module="midnightblue")
head(edges)
##              Gene1          Gene2    Weight
## 45  Zm00001d001857 Zm00001d002384 0.9401886
## 89  Zm00001d001857 Zm00001d002690 0.9675345
## 90  Zm00001d002384 Zm00001d002690 0.9185426
## 133 Zm00001d001857 Zm00001d003962 0.7178340
## 134 Zm00001d002384 Zm00001d003962 0.6534956
## 135 Zm00001d002690 Zm00001d003962 0.6840004

网络可视化

plot_gcn(
    edgelist_gcn = edges, 
    net = net, 
    color_by = "module", 
    hubs = hubs
)

网络统计：

函数net_stats()可用于计算主要网络统计数据（或属性或指数），即：connectivity、scaled connectivity、clustering coefficient、maximum adjacency ratio、density、centralization、heterogeneity、diameter、betweenness、closeness,这可能需要很长时间。

就分享到这里了

可以试试

END~

仅供粉丝老铁们参考，如有侵权或错误，请联系删除改正～

收录于合集 #生信下游分析

46个

上一篇Searchlight2：批量 RNA-seq 可视化和分析管道

喜欢此内容的人还喜欢

微信扫一扫
关注该公众号