R tips：cellphonedb细胞通讯网络图优化

生信菜鸟团

发布于 2024-07-22 19:53:31

文章被收录于专栏：生信菜鸟团

上次讲到的使用cellphonedb进行细胞通讯分析，其中的网络图的效果不是特别好，本文会就网络图进行两个优化：

（1）自身互作的环形边会绕着网络中心点向外发射状分布；

（2）仅展示一个细胞发出的细胞互作时，不需要从互作数据重新生成网络，可以有更简单的方式。

调整环形边的位置

使用公共数据集进行网络图绘制：https://github.com/elliefewings/cellphonedb_shiny/blob/master/example_cellphoneDB/pvalues.txt 这个数据结果看起来应该是cellphonedb V3的结果。

参考上次的推文，略作修改，结果如下：

library(tidyverse)
library(igraph)

# 解析cellphonedb输出文件：means文件
network_dat <- "pvalues.txt" %>% 
  read_tsv() %>% 
  dplyr::select(interacting_pair, where(is.numeric)) %>%
  pivot_longer(
    -interacting_pair,
    names_to = 'cell_pair',
    values_to = 'p'
  ) %>% 
  dplyr::filter(p < 0.05) %>% 
  group_by(cell_pair) %>% 
  summarise(counts = n()) %>% 
  mutate(SOURCE = str_split(cell_pair, "\\|") %>% map_chr(1),
         TARGET = str_split(cell_pair, "\\|") %>% map_chr(2)) %>% 
  dplyr::select(SOURCE, TARGET, counts) %>% 
  dplyr::filter(counts != 0) 



# 整体网络图

all_nodes <- unique(c(network_dat$SOURCE, network_dat$TARGET))
# [1] "Myeloid"   "NKcells_0" "NKcells_1" "Tcells"   


# 定义颜色
# 选取Dark2配色的四种颜色，分配到四种细胞上
cell_color <- RColorBrewer::brewer.pal(length(all_nodes), "Dark2")
names(cell_color) <- all_nodes


net <- graph_from_data_frame(network_dat)

E(net)$width <- E(net)$counts  %>% {((. - min(.))/diff(range(.))) * 3 + 1}
# E(net)$label <- E(net)$counts

color_match_table     <- cell_color[names(cell_color) %in% all_nodes] 
E(net)$color <- attr(E(net), 'vnames' ) %>% str_remove("\\|.*$")  %>% { color_match_table[.]}
V(net)$color <- color_match_table[names(V(net))] 


# plot

karate_groups <- cluster_optimal(net)
coords        <- layout_in_circle(net, order = order(membership(karate_groups)))  # 设置网络布局

plot(
  net,
  edge.curved      =  0.3, # 线条的弯曲程度
  edge.arrow.size  =  0.3, # 线条末端箭头大小
  margin           =  rep(0.1, 4), # 调整网络图的margin，如果字体显示不全，可以适当调大此值
  layout           =  coords, # 端点坐标
  main             =  "cellphonedb"
)

绘制的网络图如下：

rtips-network-1

可以看到这个网络图的环形边都是朝向右侧的，本文需要解决这个问题。

题外话，先聊一聊这里的layout，也就是上面传入plot的coords变量。

igraph里面有很多的布局，查看layout的帮助文档就可以知道：?igraph::layout。

比如，我们可以绘制grid布局和tree形布局：

# 一行两列绘图
opar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))

# grid
plot(
  net,
  edge.curved      =  0.3, # 线条的弯曲程度
  edge.arrow.size  =  0.3, # 线条末端箭头大小
  margin           =  rep(0.1, 4), # 调整网络图的margin，如果字体显示不全，可以适当调大此值
  layout           =  layout_on_grid(net),  # 设置网络布局
  main             =  "grid"
)

# tree
plot(
  net,
  edge.curved      =  0.3, # 线条的弯曲程度
  edge.arrow.size  =  0.3, # 线条末端箭头大小
  margin           =  rep(0.1, 4), # 调整网络图的margin，如果字体显示不全，可以适当调大此值
  layout           =  layout_as_tree (net),  # 设置网络布局
  main             =  "tree"
)

结果如下：

rtips-network-2

而且我们也可以自己手动计算坐标传给igraph，比如绘制一个正三角形的layout，并分别把myeloid和Tcells放置到中心点：

# 先查看端点的名字
attr(V(net), 'names')
# [1] "Myeloid"   "NKcells_0" "NKcells_1" "Tcells"  

# myeloid在原点 -----------
coords_tri <- matrix(c(
    c(0, 0), # 第一个点是myeloid，放置到原点
    c(0, (sqrt(3) - 1/sqrt(3))),
    c(-1, - 1/sqrt(3)),
    c(1, - 1/sqrt(3))
  ),
  ncol = 2, byrow = TRUE
)

plot(
  net,
  edge.curved      =  0.3, # 线条的弯曲程度
  edge.arrow.size  =  0.3, # 线条末端箭头大小
  margin           =  rep(0.1, 4), # 调整网络图的margin，如果字体显示不全，可以适当调大此值
  layout           =  coords_tri,  # 设置网络布局
  main             =  "triangle - 1"
)

# Tcells在原点 -----------
coords_tri <- matrix(c(
  c(0, (sqrt(3) - 1/sqrt(3))),
  c(-1, - 1/sqrt(3)),
  c(1, - 1/sqrt(3)),
  c(0, 0) # Tcells在原点
),
ncol = 2, byrow = TRUE
)

plot(
  net,
  edge.curved      =  0.3, # 线条的弯曲程度
  edge.arrow.size  =  0.3, # 线条末端箭头大小
  margin           =  rep(0.1, 4), # 调整网络图的margin，如果字体显示不全，可以适当调大此值
  layout           =  coords_tri,  # 设置网络布局
  main             =  "triangle - 2"
)

# 恢复一行一列绘图
par(opar)

结果如下：

rtips-network-3

好的，关于layout的题外话聊完了。

接着解决环形边的问题，可以看到上面图的自交互的环形边都是统一朝向右侧的，我们现在需要他们绕着中心点向外发散布局。

其原理如下图，其实就是根据每个端点的坐标，求其y/x的反正切函数值，然后将环形边的loop.angle参数设置为这个角度即可。

有两个细节：

1 igraph设置的loop.angle是顺时针计算的，所以需要将将这个角度值取相反数；

2 arctan的值域是-1/2pi到1/2pi，在第一四象限，这个角度值都是可以的，如果是左侧的二三象限，需要加上pi，反转一下方向。

rtips-network-4

具体实现如下：

###  调整节点位置的线条角度

# 添加一下互作数量
E(net)$label <- E(net)$counts

# 先找到环形边的索引，特征是起点和终点是一致的
edge.start <- igraph::ends(net, es=igraph::E(net), names=FALSE)
loop.edge.idx <- which(edge.start[,2] == edge.start[,1])

# 计算反正切函数值，并还相反数和加pi处理
angles <- atan(coords[igraph::V(net),2]/coords[igraph::V(net),1])
loop.angle <- ifelse(coords[igraph::V(net), 1]>0, -angles, pi-angles)

# 只调整环形边的旋转角度
igraph::E(net)$loop.angle <- 0 # loop.angle的默认值就是0，但是基于稳健性考虑，此句不可省略
igraph::E(net)$loop.angle[loop.edge.idx] <- loop.angle[edge.start[loop.edge.idx, 1]    ]

plot(
  net,
  edge.curved      = 0.3,
  edge.arrow.size  = 0.3,
  margin           = rep(0.3, 4), # 有遮挡，适当调大margin
  layout           = coords,
  main             = "cellphonedb - adjust loop edge"
)

结果如下：

rtips-network-5

简化绘制拆分的环形网络图

由于网络图的绘制元素中，只需要将特定的边或者文字颜色置空就可以将其取消绘制，所以拆分的环形网络图并不需要从头绘制，只需要取消特定元素的绘制即可。

all_nodes <- attr(V(net), "names")
all_edges <- attr(E(net), "vnames")

# 两行两列绘图
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar = c(1, 0))

all_nodes %>%
  walk(function(node){
    net2 <- net

    # 当前node为起点的边
    edge_idx <- str_starts(all_edges, fixed(node))

    # 只保留当前节点发出的线和文字
    E(net2)$color[! edge_idx] <- NA
    E(net2)$label[! edge_idx] <- NA


    plot(
      net2,
      edge.curved      = 0.3,
      edge.arrow.size  = 0.3,
      margin           = rep(0.3, 4),
      layout           = coords,
      main             = node
    )
  })

# 绘图一行一列绘图
par(opar)

结果如下：