search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
127 stars 32 forks source link

MEGENA网络分析 #1857

Closed ixxmu closed 2 years ago

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/n-CdXu1397Pqps8y_MNDxQ

github-actions[bot] commented 2 years ago

MEGENA网络分析 by 生信技能树


给学徒们收集整理了几套带GitHub源代码的文献图表合辑，让优秀者一点一滴拆解开来分享给大家。（全部的代码复制粘贴即可运行，欢迎尝试以及批评指正）

下面是去年实习生的分享
 author: "ylchen" 

WGCNA大家应该是已经很熟悉了，我们已经在生信技能树多个教程分享WGCNA的实战细节，见：

  • 一文学会WGCNA分析

  • 一文看懂WGCNA 分析(2019更新版) （点击阅读原文即可拿到测序数据）

  • 通过WGCNA作者的测试数据来学习

  • 重复一篇WGCNA分析的文章（代码版）

  • 重复一篇WGCNA分析的文章（解读版）（逆向收费读文献2019-19）

  • 关键问题答疑：WGCNA的输入矩阵到底是什么格式

  • WGCNA-流程及原理细节直播互动授课（今晚八点）

其他：

但是如果你想要一些新意，那不如看看MEGENA(多尺度嵌入基因共表达网络分析)？

image-20220223204244635

这个网络是由来自伊坎西奈山医学院的科学家们开发的，并于2015年发表《PLoS Computational Biology》中。目前已应用于多种疾病如癌症、肥胖，阿尔茨海默氏症。测试数据来自癌症基因组图谱(TCGA)，并且MEGENA在乳腺癌和肺癌中鉴定了新的监管目标，优于其它共表达分析方法。

其实虽然两者也是根据基因-基因相关性分成模块，求与疾病最相关的模块出来，进而筛选关键基因。不过两者肯定是用不同的统计学方法，有兴趣的小伙伴可以比较以下。

Multiscale Embedded Gene Co-expression Network Analysis (plos.org)

MEGENA.pdf (r-project.org)

MEGENA_pipeline (r-project.org)

下面介绍下MEGENA包的主要用法

0.前言

MEGENA算法的一个关键部分是构造平面滤波网络(PFN)。包括三个主要步骤：1）根据基因表达谱进行相关性分析，筛选出关键基因对构建平面滤波网络(PFN)；2）从PFN开始，对模块间进行多尺度聚类分析，得到层次聚类结果；3）采用多尺度中心分析（MHA）来识别高度复杂的网络聚类特征（CTA），探讨模块与临床结果的相关性。

图1显示了MEGENA的总体分析流程。

Untitled

1. R包安装

rm(list = ls()) # rm R working space
# library(devtools);
# install_github("songw01/MEGENA");
# BiocManager::install('MEGENA')
library(MEGENA)

2.数据载入

  • 数据集：datExpr
  • A portion of TCGA breast cancer data to test run MEGENA（a gene expression matrix.）
# load toy example data
data(Sample_Expression) 
dim(datExpr) # TCGA breast cancer data,a gene expression matrix
datExpr[1:4,1:4]
> dim(datExpr)
[1] 300 844
> datExpr[1:4,1:4]
         TCGA-E9-A1RC-01A-11R-A157-07 TCGA-BH-A0B1-01A-12R-A056-07
C1QB                         7.799198                     9.491122
CD3E                         2.703622                     6.570632
CD3D                         1.121214                     4.802647
TMEM176A                     6.608023                     8.948846
         TCGA-B6-A0RG-01A-11R-A056-07 TCGA-D8-A1JH-01A-11R-A13Q-07
C1QB                         9.959445                    10.648426
CD3E                         5.812081                     8.074163
CD3D                         4.104806                     6.152895
TMEM176A                    10.054326                     8.760917

3.计算基因表达谱之间的相似性

  • 函数：calculate.correlation
  • correlation analysis with FDR calculation

原理：用错误发现率(fdr)筛选基因相似性，以尽量减少错误阳性的影响。使用平面滤波网络(PFN)计算有显著相关性的基因对  ,在使用不同相似度度量的情况下，可以独立地将结果格式化为3列的数据帧，列名为c(row, col, weight)，并确保权重列在0到1之间

# input parameters参数
n.cores <- 2; # 内核/线程数
doPar <-TRUE; # 是否需要并行
method = "pearson" # 计算相关性的方法pearson or spearman. 
FDR.cutoff = 0.05 # FDR阈值，定义相关性
module.pval = 0.05 # 模块显著性，默认0.05. 
hub.pval = 0.05 # 随机四面体网络连接度p值
cor.perm = 10; # number of permutations for calculating FDRs for all correlation pairs. 
hub.perm = 100; # number of permutations for calculating connectivity significance p-value.
# correlation analysis with FDR calculation
ijw <- calculate.correlation(datExpr,
                             doPerm = cor.perm,
                             output.corTable = FALSE, # 输出表格
                             output.permFDR = FALSE)
head(ijw)
## i = 1
## i = 2
## i = 3
## i = 4
## i = 5
## i = 6
## i = 7
## i = 8
## i = 9
## i = 10
> head(ijw)
          row  col       rho
C1QB.127 C1QB C1QC 0.9793833
CD3E.6   CD3E  CD2 0.9793832
CD3D.5   CD3D  CD2 0.9745410
C1QB.3   C1QB C1QA 0.9734056
CD3E     CD3E CD3D 0.9730448
C1QA.123 C1QA C1QC 0.9628771

4. 构建PFN

  • 函数：calculate.PFN
  • main function to calculate PFN a ranked list of edge pairs
##### calculate PFN
el <- calculate.PFN(ijw[,1:3],
                    doPar = doPar,# 并行
                    num.cores = n.cores, # 核数
                    keep.track = FALSE) # 创建pfg_el.RData
## ####### PFN Calculation commences ########
## [1] "343 out of 894 has been included."
## [1] "471 out of 894 has been included."
## [1] "556 out of 894 has been included." 
## [1] "PFG is complete."
g <- graph.data.frame(el,directed = FALSE) #转化为图数据结构
# directed = FALSE为无向网络

5.MCA聚类

  • 函数：do.MEGENA
  • multiscale clustering analysis (MCA) and multiscale hub analysis (MHA) pipeline

通过MCA聚类进行多尺度聚类分析。“MEGENA.output “是用于下游分析的核心输出，用于汇总和绘图

##### perform MCA clustering.
MEGENA.output <- do.MEGENA(g,
                           mod.pval = module.pval, # 模块显著性，默认0.05. 
                           hub.pval = hub.pval,# 随机四面体网络连接度p值
                           remove.unsig = TRUE, # 是否移除不显著的cluster
                           min.size = 10,# 最小cluster数
                           max.size = vcount(g)/2, # 最大cluster数
                           doPar = doPar, # 并行
                           num.cores = n.cores, # 核数
                           n.perm = hub.perm, # 计算连接度P值的排列数
                           save.output = FALSE) # 保存每一步分析

6. 汇总结果

  • 函数：MEGENA.ModuleSummary
  • Summarizes modules into a table
# annotation to be done on the downstream
annot.table=NULL
id.col = 1
symbol.col= 2
###########
summary.output <- MEGENA.ModuleSummary(MEGENA.output,
                                       mod.pvalue = module.pval,hub.pvalue = hub.pval,
                                       min.size = 10,max.size = vcount(g)/2,
                                       annot.table = annot.table,id.col = id.col,symbol.col = symbol.col,
                                       output.sig = TRUE)

这里要注意annot.table，不过看说明

 

应该就是节点名≠基因名时，可以再载入个注释表格进行转换，不过指定是哪两列

summary.output$modules

结果得到13个模块

> str(summary.output$modules)
List of 13
 $ c1_2 : chr [1:133] "C1QB" "C1QA" "SASH3" "HLA-DPA1" ...
 $ c1_3 : chr [1:121] "CD3E" "CD3D" "CD2" "CXCL11" ...
 $ c1_4 : chr [1:22] "FCRL5" "CD79A" "LOC96610" "POU2AF1" ...
 $ c1_5 : chr [1:24] "IFIT3" "IFI44L" "IFIT1" "IFI6" ...
 $ c1_6 : chr [1:26] "C1QB" "C1QA" "FCGR1A" "FCGR1C" ...
 $ c1_7 : chr [1:68] "SASH3" "PLEK" "NCKAP1L" "EVI2B" ...
 $ c1_10: chr [1:65] "CD3E" "CD3D" "CD2" "ITK" ...
 $ c1_12: chr [1:25] "SLAMF6" "LY9" "GPR174" "CD38" ...
 $ c1_13: chr [1:11] "TAP1" "HLA-B" "HLA-F" "HLA-H" ...
 $ c1_15: chr [1:16] "CD79A" "MS4A1" "CD19" "BLK" ...
 $ c1_19: chr [1:22] "SASH3" "ARHGAP9" "PTPN7" "CD48" ...
 $ c1_21: chr [1:56] "CD3E" "CD3D" "CD2" "ITK" ...
 $ c1_23: chr [1:11] "SLAMF6" "LY9" "ZNF831" "BTLA" ...

summary.output$module.table

Untitled

若提供了"annot.table"，还会出现gene symbol和id对应的结果（mapped.modules）

> print(summary.output$module.table)
      module.id module.size module.parent
c1_2       c1_2         133          c1_1
c1_3       c1_3         121          c1_1
c1_4       c1_4          22          c1_1
c1_5       c1_5          24          c1_1
c1_6       c1_6          26          c1_2
c1_7       c1_7          68          c1_2
c1_10     c1_10          65          c1_3
c1_12     c1_12          25          c1_3
c1_13     c1_13          11          c1_3
c1_15     c1_15          16          c1_4
c1_19     c1_19          22          c1_7
c1_21     c1_21          56         c1_10
c1_23     c1_23          11         c1_12
                                                                         module.hub
c1_2  SASH3(34),CD53(33),PTPRC(26),CD86(20),NCKAP1L(19),FERMT3(19),PLEK(18),CD4(16)
c1_3                                   CD2(33),CD3E(32),GBP5(20),ITK(19),SLAMF6(18)
c1_4                                                                      MS4A1(12)
c1_5                                                                      IFIT3(19)
c1_6                                                                       CD86(14)
c1_7                                CD53(32),SASH3(27),PLEK(17),NCKAP1L(17),CD4(16)
c1_10                                                              CD2(32),CD3E(29)
c1_12                                                                    SLAMF6(14)
c1_13                                                                      PSMB9(9)
c1_15                                                           MS4A1(12),FCRL1(11)
c1_19                                                           SASH3(19),PTPN7(12)
c1_21                                                              CD2(32),CD3E(29)
c1_23                                                                    SLAMF6(10)
      module.scale module.pvalue
c1_2            S4          0.00
c1_3          <NA>          0.00
c1_4            S4          0.00
c1_5            S4          0.00
c1_6            S3          0.01
c1_7          <NA>          0.00
c1_10         <NA>          0.00
c1_12           S4          0.00
c1_13           S4          0.00
c1_15           S3          0.03
c1_19           S3          0.00
c1_21           S4          0.00
c1_23           S3          0.00

7.可视化

包里还自带绘图函数

pnet.obj <- plot_module(output.summary = summary.output,
                        PFN = g,
                        subset.module = "c1_4", # 可视化的cluster
                        layout = "kamada.kawai", # 布局的算法:"kamada.kawai"or"fruchterman.reing"
                        label.hubs.only = FALSE, # 只展示关键基因的标签
                        color.code =  "grey",
                        output.plot = FALSE, # 输出图片
                        out.dir = "modulePlot", # 新建文件夹存结果
                        col.names = c("magenta","green","cyan"), # 网络的颜色
                        label.scaleFactor = 2, # 节点标签大小
                        hubLabel.col = "black", # 标签颜色
                        hubLabel.sizeProp = 1, # 调整关键节点标签与普通的比例，默认0.5
                        show.topn.hubs = Inf, # 显示的最大关键节点数量
                        show.legend = TRUE) # 注释
print(pnet.obj[[1]])
Untitled

因为这个网络是有层次结果的，所以输出来的图会呈现出父模块与子模块间的网络关系

不过做出来的图挺丑的，那还不如用老办法，把数据输入到cytoscape做个漂亮的出来

sub_cluster <- summary.output$modules$c1_4 # 模块基因
# 与前面算出权重的el数据集取交集
df1 = el[el[,1] %in% sub_cluster,]
df2 = df1[df1[,2] %in% sub_cluster,]
# df2就可以直接输入cytoscape了，有节点，有边的关系
> head(df2)
         row      col    weight
22     FCRL5  POU2AF1 0.9517072
62     CD79A  POU2AF1 0.9370254
97     FCRL5 LOC96610 0.9295689
105    FCRL5    CD79A 0.9275319
129 LOC96610 TNFRSF17 0.9230158
138    FCRL5 TNFRSF17 0.9218706

8.关键还是得多看文献

不过具体MEGENA网络的用法，如何结合WGCNA网络挖掘关键生物信息呢？关键还是得多看文献哦！

下面近五年用了MEGENA网络方法的高分文章，可以从文献中找到答案哦

Quiescent stem cell marker genes in glioma gene networks are sufficient to distinguish between normal and glioblastoma (GBM) samples | Scientific Reports (nature.com)

Network models of primary melanoma microenvironments identify key melanoma regulators underlying prognosis | Nature Communications

Multiscale network analysis reveals molecular mechanisms and key regulators of the tumor microenvironment in gastric cancer - Song - 2020 - International Journal of Cancer - Wiley Online Library

Network models of prostate cancer immune microenvironments identify ROMO1 as heterogeneity and prognostic marker | Scientific Reports (nature.com)

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 jmzeng1314@163.com

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你