在前幾篇里我用cytoscape軟件構(gòu)建了網(wǎng)絡(luò)圖（如何將WGCNA導(dǎo)出的基因共表達(dá)網(wǎng)絡(luò)進(jìn)行可視化），見(jiàn)下圖。但是看上去好像什么信息也看不出來(lái)，因?yàn)榛蚝突蛑g的connectivity實(shí)在是太多了，全都糊在一起，看不出有哪些基因位于網(wǎng)絡(luò)的中心，或者說(shuō)看不出哪些基因和其他基因的聯(lián)系很多。
當(dāng)然你可以把cytoscape里的每一個(gè)基因的degree.layout導(dǎo)出來(lái)，里面有每一個(gè)基因有多少條edge與其相連，然后做一個(gè)bar圖（用R畫Circular barplot圖）,但是對(duì)于圖里的數(shù)字（比如60和20）并沒(méi)有非常直觀的展示。偶然看到簡(jiǎn)書里一篇文章寫的非常好（R數(shù)據(jù)可視化18:弧形圖）,做出的圖形是這樣的：
這樣的圖可以清晰的顯示出哪些基因之間的聯(lián)系非常多，哪些基因形成一個(gè)“小社區(qū)”。圖片來(lái)自：http://www.itdecent.cn/p/a236227503e4
（一）弧形圖需要的數(shù)據(jù)
弧形圖需要的數(shù)據(jù)是從cytoscape里導(dǎo)出的，你也可以自己制作（如果數(shù)據(jù)不多的話）：
在cytoscape右下角處選擇“edge table”，然后點(diǎn)擊“export”，就會(huì)得到一個(gè)excel表
導(dǎo)出的表長(zhǎng)這樣：
我們需要的是這個(gè)表的“name”列。你需要在excel里把“ (interacts with) ”替換成空格（注意括號(hào)兩邊也有空格），這一步在excel里相當(dāng)簡(jiǎn)單，就不贅述了。替換完成后的表長(zhǎng)這樣：
現(xiàn)在基因和基因之間只有一個(gè)空格，一會(huì)兒要在R里根據(jù)這個(gè)空格將這一列分成兩列
（二）數(shù)據(jù)的預(yù)處理
> x = read.csv("turquiose_weight_morethan0.35_edge_interaction_edited.csv",header = T,sep = ",")
#先從原始矩陣?yán)锾崛⌒畔ⅲ?gòu)建需要的數(shù)據(jù)框
> data = as.data.frame(x$name)
> colnames(data) = "V"
#把這列按照空格分成2列
> library(stringr)
> data <- str_split_fixed(data$V," ", 2)
> colnames(data) = c("from","to")
#給表格加一列，這一列是注明兩個(gè)基因間是否有作用，如果有作用，標(biāo)注1，沒(méi)有的話標(biāo)注0
#這里每一行都是2個(gè)基因，所以我們要加一列“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
現(xiàn)在得到的表格就是我們需要的。
（三）計(jì)算基因的interact次數(shù)（edge條數(shù)）
> library(tidyverse)
> c(as.character(data$from), as.character(data$to)) %>%
as_tibble() %>%
group_by(value) %>% #按照value進(jìn)行分組
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
（四）給基因鑒定“小社區(qū)”
> library(igraph)
> mygraph <- graph_from_data_frame(data, vertices = co, directed = FALSE)
# 鑒定社區(qū)
> com <- walktrap.community(mygraph)
#重新整理數(shù)據(jù)
> 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")) #調(diào)色板請(qǐng)看：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))
