厚缊

诹图——相关性网络

厚缊 / 2020-01-14


时间过得太快,眼瞅着又过了半个月,ggcor 0.9.2 版本终于完工了,还剩一个小时坐车回家,写篇文章结束2019。需要事先说明,ggcor仅仅打通了和其它R包和软件的数据联系,自身没有提供任何网络可视化的工具。

安装

尚处于测试阶段,不推荐在正式场合使用。

devtools::install_github("houyunhuang/ggcor")

相关系数矩阵转化为网络

ggcor提供了cor_network()as_graph_tbl()函数来进行网络数据的转化,支持psych::corr.test()Hmisc::rcorr()ggcor::correlate()ggcor::fast_correlate()函数相关性分析结果直接转成graph_tblggraph网络可视化的数据结构),并且能同时完成网络边和节点的过滤,让相关性网络构建和可视化大大简化。

cor_network()函数

cor_network()函数是ggcor数据过滤和网络转化的基础。corrp.value是相关系数及检验p值矩阵,其它重要参数是r.thresr.absoluter.thres,当r.absolute为TRUE,表示按照相关系数的绝对值进行过滤,r.thresp.thres是相关系数、统计显著性过滤的门槛值。cor_network()会自动去除AB-BA形式的重复值,前提是相关系数矩阵都是有行名和列名的。

library(ggcor)
args(cor_network)
## function (corr, p.value = NULL, row.names = NULL, col.names = NULL, 
##     simplify = TRUE, r.thres = 0.6, r.absolute = TRUE, p.thres = 0.05) 
## NULL
cor(mtcars) %>% cor_network()
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 3
##    from    to      r
##   <int> <int>  <dbl>
## 1     1    11 -0.852
## 2     2    11 -0.848
## 3     1     2  0.902
## # … with 23 more rows

as_graph_tbl()函数

fortify_cor(mtcars, type = "upper") %>% 
  as_tbl_graph() 
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 5
##    from    to      r .row.id .col.id
##   <int> <int>  <dbl>   <int>   <int>
## 1     1    11 -0.852      11       2
## 2     2    11 -0.848      11       3
## 3     1     2  0.902      10       3
## # … with 23 more rows
correlate(mtcars) %>% 
  as_tbl_graph() 
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 3
##    from    to      r
##   <int> <int>  <dbl>
## 1     1    11 -0.852
## 2     2    11 -0.848
## 3     1     2  0.902
## # … with 23 more rows
fast_correlate(mtcars) %>% 
  as_tbl_graph()
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 4
##    from    to      r  p.value
##   <int> <int>  <dbl>    <dbl>
## 1     1    11 -0.852 6.11e-10
## 2     2    11 -0.848 9.38e-10
## 3     1     2  0.902 1.80e-12
## # … with 23 more rows
Hmisc::rcorr(data.matrix(mtcars)) %>% 
  as_tbl_graph()
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 4
##    from    to      r  p.value
##   <int> <int>  <dbl>    <dbl>
## 1     1    11 -0.852 6.11e-10
## 2     2    11 -0.848 9.38e-10
## 3     1     2  0.902 1.80e-12
## # … with 23 more rows
psych::corr.test(mtcars) %>% 
  as_tbl_graph()
## # A tbl_graph: 11 nodes and 26 edges
## #
## # An undirected simple graph with 1 component
## #
## # Node Data: 11 x 1 (active)
##   name 
##   <chr>
## 1 cyl  
## 2 disp 
## 3 hp   
## 4 drat 
## 5 wt   
## 6 qsec 
## # … with 5 more rows
## #
## # Edge Data: 26 x 4
##    from    to      r  p.value
##   <int> <int>  <dbl>    <dbl>
## 1     1    11 -0.852 3.18e- 8
## 2     2    11 -0.848 4.78e- 8
## 3     1     2  0.902 9.92e-11
## # … with 23 more rows

网络可视化

有了以上两个函数,基本能打通与所有R网络可视化包之间的联系。之所以如此肯定,是因为R中网络可视化算法的两大基础包igraphsna本身是相互联系的,而graph_tbl可以很方便的转成igraph

library(vegan)
library(tidygraph)
library(ggraph)
library(igraph)
data("varespec")
correlate(varespec) %>% 
  as_tbl_graph() %>% 
  ggraph("circle") + 
  geom_edge_fan(aes(edge_width = r), edge_colour = "grey60",
                show.legend = FALSE) +
  geom_node_point(aes(colour = name), size = 5, show.legend = FALSE) +
  scale_edge_width_continuous(range = c(0.1, 2)) +
  coord_fixed() +
  theme_graph()

correlate(varespec) %>% 
  as_tbl_graph() %>% 
  as.igraph() %>% 
  plot(layout = layout.circle)

处理效率

对于生物信息领域,进行相关性分析的数据往往非常大,一般的R包超过2000个变量基本很难运行。目前来看,WGCNA是领域计算速度最快的,处理5000个变量压力不大。ggcor本身不依赖WGCNA包,要使用以下功能,请先安装。

我的电脑处理50个样本,5000个变量的运行时间是4秒左右,应该算是可以接受的时间了。

m <- matrix(rnorm(50 * 5000), nrow = 50)
rownames(m) <- paste0("row", 1:50)
colnames(m) <- paste0("col", 1:5000)
corr <- fast_correlate(m)
microbenchmark::microbenchmark(as_tbl_graph(corr), times = 5)
## Unit: seconds
##                expr      min       lq     mean   median       uq      max
##  as_tbl_graph(corr) 2.970127 3.223403 3.378703 3.230509 3.585272 3.884205
##  neval
##      5

数据导出

网络图确实使用gephi,cytoscape更简单,我们可以用ggcor做数据预处理,导出后用专门的网络可视化工具绘图。graph_tbl能直接转成tibble进行导出。

导出边数据

edges <- Hmisc::rcorr(data.matrix(varespec)) %>% 
  as_tbl_graph() %>% 
  activate("edges") %>% 
  as_tibble()
head(edges)
## # A tibble: 6 x 4
##    from    to     r     p.value
##   <int> <int> <dbl>       <dbl>
## 1     1    25 0.676 0.000286   
## 2     2    26 0.611 0.00151    
## 3     3    25 0.610 0.00156    
## 4     4    25 0.821 0.000000898
## 5     5    25 0.660 0.000444   
## 6     4     5 0.851 0.000000139

导出节点数据

nodes <- Hmisc::rcorr(data.matrix(varespec)) %>% 
  as_tbl_graph() %>% 
  activate("nodes") %>% 
  as_tibble()
head(nodes)
## # A tibble: 6 x 1
##   name    
##   <chr>   
## 1 Vaccmyrt
## 2 Vaccviti
## 3 Descflex
## 4 Betupube
## 5 Dicrpoly
## 6 Hylosple