用R画弧形图（Arc Diagram）

在前几篇里我用cytoscape软件构建了网络图（如何将WGCNA导出的基因共表达网络进行可视化），见下图。但是看上去好像什么信息也看不出来，因为基因和基因之间的connectivity实在是太多了，全都糊在一起，看不出有哪些基因位于网络的中心，或者说看不出哪些基因和其他基因的联系很多。

当然你可以把cytoscape里的每一个基因的degree.layout导出来，里面有每一个基因有多少条edge与其相连，然后做一个bar图（用R画Circular barplot图）,但是对于图里的数字（比如60和20）并没有非常直观的展示。偶然看到简书里一篇文章写的非常好（R数据可视化18:弧形图）,做出的图形是这样的：

这样的图可以清晰的显示出哪些基因之间的联系非常多，哪些基因形成一个“小社区”。图片来自：https://www.jianshu.com/p/a236227503e4

（一）弧形图需要的数据

弧形图需要的数据是从cytoscape里导出的，你也可以自己制作（如果数据不多的话）：

在cytoscape右下角处选择“edge table”，然后点击“export”，就会得到一个excel表

导出的表长这样：

我们需要的是这个表的“name”列。你需要在excel里把“ (interacts with) ”替换成空格（注意括号两边也有空格），这一步在excel里相当简单，就不赘述了。替换完成后的表长这样：

现在基因和基因之间只有一个空格，一会儿要在R里根据这个空格将这一列分成两列

（二）数据的预处理

> x = read.csv("turquiose_weight_morethan0.35_edge_interaction_edited.csv",header = T,sep = ",")
#先从原始矩阵里提取信息，并构建需要的数据框
> data = as.data.frame(x$name)
> colnames(data) = "V"
#把这列按照空格分成2列
> library(stringr)
> data <- str_split_fixed(data$V," ", 2)
> colnames(data) = c("from","to")
#给表格加一列，这一列是注明两个基因间是否有作用，如果有作用，标注1，没有的话标注0
#这里每一行都是2个基因，所以我们要加一列“1”
> to_add = as.data.frame(rep(1,1245)) #每一行都是1，一共是1245行
> colnames(to_add) = "value"
> data = cbind(data,to_add)
> head(data)
    from    to value
1  Spry1 Fbln1     1
2 Pdgfrl Fbln1     1
3    Dcn  Rgl1     1
4  Spry1 Prelp     1
5  Spry1    Cp     1
6  Spry1  Tgm2     1

现在得到的表格就是我们需要的。

（三）计算基因的interact次数（edge条数）

> library(tidyverse)
> c(as.character(data$from), as.character(data$to)) %>%
  as_tibble() %>% 
  group_by(value) %>% #按照value进行分组
  summarize(n=n()) -> co
#修改列名
> colnames(co) <- c("name", "n")
#排序，按照降序排列
> co <- co[order(co$n,decreasing = T),]
#分组(for color)
> A = sum(co$n >=60)
> B = sum(co$n >=50) - A 
> C = sum(co$n >=40) - A - B
> D = sum(co$n >=30) - (A+B+C)
> E = sum(co$n >=20) - (A+B+C+D)
> f = sum(co$n >=10) - (A+B+C+D+E)
> G = sum(co$n < 10)
#把分组信息加到表里
> co <- data.frame(
  co,
  group=c(rep('A', A), rep('B', B), rep('C', C), rep('D', D),rep('E',E),rep('F',f),rep('G',G))
)
> head(co) 
   name  n group
1 Spry1 68     A
2  C1ra 67     A
3 Hmox1 67     A
4 Timp3 55     B
5  Rgl1 54     B
6 Prelp 53     B

（四）给基因鉴定“小社区”

> library(igraph)
> mygraph <- graph_from_data_frame(data, vertices = co, directed = FALSE)
# 鉴定社区
> com <- walktrap.community(mygraph)
#重新整理数据
> coline <- co %>% 
  mutate(grp = com$membership) %>%
  arrange(grp) %>%
  mutate(name=factor(name, name))
> mygraph2 <- graph_from_data_frame(data, vertices = coline, directed = FALSE)

（五）作图

#作图
> library(RColorBrewer)
> library(ggraph)
> mycolor<-colorRampPalette(brewer.pal(7, "Set3"))  #调色板请看：https://www.datanovia.com/en/blog/the-a-z-of-rcolorbrewer-palette/
> ggraph(mygraph2, layout="linear") + 
  geom_edge_arc(edge_colour="black", edge_alpha=0.2, edge_width=0.3, fold=TRUE) +
  geom_node_point(aes(size=n, color=as.factor(grp), fill=co$group), alpha=0.6) +
  scale_size_continuous(range=c(1,10)) +
  scale_color_manual(values=mycolor(50)) +
  geom_node_text(aes(label=name), angle=90, hjust=1, nudge_y = -1.1, size=2.8) +
  theme_void() +
  theme(
    legend.position="none",
    plot.margin=unit(c(0,0,0.5,0), "null"),
    panel.spacing=unit(c(0,0,2,0), "null") ) +
  expand_limits(x = c(-3, 1.2), y = c(-14, 1.2))