之前我們微生物生態網絡基於igraph進行出圖，這是基礎包寫的出圖函數，近年來，隨着ggplot的普及，使用ggplot出圖似乎成了部分人的基本功，plot函數為代表的基本包已經很少有人去學習了。
為了減少學習成本，同時制作多樣化程度和個性化程度更高的網絡圖，使用ggplot出圖成了必然的選擇，近年來很多人基於ggplot的版本的網絡圖進行了許多嘗試，但是就微生物網絡而言，這一方面使用基於ggplot為基礎網絡的人很少。
一方面因為轉移成本問題，一方面因為展示樣式問題。 慢慢的工作量積累夠了：
16年gplot.layout的出現，使得ggplot網絡可以擴展19中可視化方式。
ggplot出圖數據為矩陣，依托於強大的數據框處理函數aplyr包，我們可以對網絡進行個性化程度極高的設置，包括標簽，圖例，顏色，形狀，大小。
這是google上的一個嘗試
通過ggplot做出基本網絡圖形
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(network))
suppressPackageStartupMessages(library(sna))
suppressPackageStartupMessages(library(ergm))
library(network)
library(ggplot2)
library(sna)
library(ergm)
plotg <- function(net, value = NULL) {
m <- as.matrix.network.adjacency(net) # get sociomatrix
# get coordinates from Fruchterman and Reingold's force-directed placement
# algorithm.
plotcord <- data.frame(gplot.layout.fruchtermanreingold(m, NULL))
# or get it them from Kamada-Kawai's algorithm: plotcord <-
# data.frame(gplot.layout.kamadakawai(m, NULL))
colnames(plotcord) = c("X1", "X2")
edglist <- as.matrix.network.edgelist(net)
edges <- data.frame(plotcord[edglist[, 1], ], plotcord[edglist[, 2], ])
plotcord$elements <- as.factor(get.vertex.attribute(net, "elements"))
colnames(edges) <- c("X1", "Y1", "X2", "Y2")
edges$midX <- (edges$X1 + edges$X2)/2
edges$midY <- (edges$Y1 + edges$Y2)/2
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2),
data = edges, size = 0.5, colour = "grey") + geom_point(aes(X1, X2,
colour = elements), data = plotcord) + scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) + theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) + theme(panel.background = element_rect(fill = "white",
colour = NA)) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
return(print(pnet))
}
g <- network(50, directed = FALSE, density = 0.03)
classes <- rbinom(50, 1, 0.5) + rbinom(50, 1, 0.5) + rbinom(50, 1, 0.5)
set.vertex.attribute(g, "elements", classes)
g
plotg(g)
嘗試：基於ggplot出圖儲存於list后批量拼圖
ggplot拼圖函數有很多，但是這里我們批量出圖儲存於list中，這里使用下面這個函數做拼圖。 google上有人寫的：
multiplot <- function(..., plotlist=NULL, cols) {
require(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# Make the panel
plotCols = cols # Number of columns of plots
plotRows = ceiling(numPlots/plotCols) # Number of rows needed, calculated from # of cols
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
vplayout <- function(x, y)
viewport(layout.pos.row = x, layout.pos.col = y)
# Make each plot, in the correct location
for (i in 1:numPlots) {
curRow = ceiling(i/plotCols)
curCol = (i-1) %% plotCols + 1
print(plots[[i]], vp = vplayout(curRow, curCol ))
}
}
基於微生物大量的OTU，我嘗試了18中layout
大家使用cor.test計算得到的相關矩陣即可作為輸入
# 確定物種間存在相互作用關系的閾值，將相關性R矩陣內不符合的數據轉換為0
occor.r[occor.p>p.threshold|abs(occor.r)<r.threshold] = 0
library(network)
library(ggplot2)
library(sna)
library(ergm)
library(igraph)
g <- network(occor.r, directed=FALSE,vertex.attrnames=T)
summary(g)
net = g
m <- as.matrix.network.adjacency(net) # get sociomatrix
plotcord = list()
plotcord[[1]] <- data.frame(gplot.layout.fruchtermanreingold(m, NULL));names(plotcord[[1]]) = "fruchtermanreingold"
plotcord[[2]] <- data.frame(gplot.layout.kamadakawai(m, NULL));names(plotcord[[2]]) = "kamadakawai"
plotcord[[ 3]] <- data.frame(gplot.layout.adj(m, NULL));names(plotcord[[3]]) = "adj"
plotcord[[ 4]] <- data.frame(gplot.layout.circle(m, NULL));names(plotcord[[4]]) = "circle"
plotcord[[ 5]] <- data.frame(gplot.layout.circrand(m, NULL));names(plotcord[[5]]) = "circrand"
plotcord[[ 6]] <- data.frame(gplot.layout.eigen(m, NULL));names(plotcord[[6]]) = "eigen"
plotcord[[ 7]] <- data.frame(gplot.layout.fruchtermanreingold(m, NULL));names(plotcord[[7]]) = "fruchtermanreingold"
plotcord[[ 8]] <- data.frame(gplot.layout.geodist(m, NULL));names(plotcord[[8]]) = "geodist"
plotcord[[ 9]] <- data.frame(gplot.layout.hall(m, NULL));names(plotcord[[9]]) = "hall"
plotcord[[ 10]]<- data.frame(gplot.layout.kamadakawai(m, NULL));names(plotcord[[10]]) = "kamadakawai"
plotcord[[11]] <- data.frame(gplot.layout.mds(m, NULL));names(plotcord[[11]]) = "mds"
plotcord[[12]] <- data.frame(gplot.layout.princoord(m, NULL));names(plotcord[[12]]) = "princoord"
plotcord[[13 ]] <- data.frame(gplot.layout.random(m, NULL));names(plotcord[[13]]) = "random"
plotcord[[14 ]] <- data.frame(gplot.layout.rmds(m, NULL));names(plotcord[[14]]) = "rmds"
plotcord[[15 ]] <- data.frame(gplot.layout.segeo(m, NULL));names(plotcord[[15]]) = "segeo"
plotcord[[16 ]] <- data.frame(gplot.layout.seham(m, NULL));names(plotcord[[16]]) = "seham"
plotcord[[17 ]] <- data.frame(gplot.layout.spring(m, NULL));names(plotcord[[17]]) = "spring"
plotcord[[18 ]] <- data.frame(gplot.layout.springrepulse(m, NULL));names(plotcord[[18]]) = "springrepulse"
#計算花費很長時間，所以不計算了
# plotcord[[19 ]] <- data.frame(gplot.layout.target(m, NULL))
plots = list()
ii = 1
for (ii in 1:18) {
plotcor = plotcord[[ii]]
colnames(plotcor) = c("X1", "X2")
head(plotcor)
plotcor$elements <- colnames(occor.r)
edglist <- as.matrix.network.edgelist(net)
edglist = as.data.frame(edglist)
# aaaa = as.matrix.network(net)
# 構建igraph對象構建鄰接矩陣
igraph <- graph_from_adjacency_matrix(occor.r,mode="undirected",weighted=TRUE,diag=FALSE)
E(igraph)$weight
edglist$weight = E(igraph)$weight
edges <- data.frame(plotcor[edglist[, 1], ], plotcor[edglist[, 2], ])
head(edges)
edges$weight = E(igraph)$weight
##這里將邊權重根據正負分為兩類
aaa = rep("a",length(edges$weight))
for (i in 1:length(edges$weight)) {
if (edges$weight[i]> 0) {
aaa[i] = "+"
}
if (edges$weight[i]< 0) {
aaa[i] = "-"
}
}
#添加到edges中
edges$wei_label = aaa
colnames(edges) <- c("X1", "Y1","OTU_1", "X2", "Y2","OTU_2","weight","wei_label")
edges$midX <- (edges$X1 + edges$X2)/2
edges$midY <- (edges$Y1 + edges$Y2)/2
head(edges)
# library(ggrepel)
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,colour = edges$wei_label),
data = edges, size = 0.5) +
geom_point(aes(X1, X2), size=3, pch = 21, data = plotcor, fill = "#8DD3C7") + scale_colour_brewer(palette = "Set1") +
scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
labs( title = names(plotcord[[ii]])[1])+
# geom_text(aes(X1, X2,label=elements),size=4, data = plotcor)+
# discard default grid + titles in ggplot2
theme(panel.background = element_blank()) +
theme(legend.position = "none") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
theme(legend.background = element_rect(colour = NA)) +
theme(panel.background = element_rect(fill = "white", colour = NA)) +
theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet
plots[[ii]] = pnet
}
pdf(file = "cs.pdf",width = 12,height = 18)
multiplot(plotlist=plots,cols=3)
dev.off()
到這里我們就可以將微生物生態網絡移植到ggolot中
這里我選擇合適和layout布局方式，使用google網上構造的拼圖工具，結合ggolot的微生物生態網絡圖 這里我做了細菌和真菌的五個處理的網絡，隱去標簽。
備注：代碼略長，此處略去；
歡迎加入微生信生物討論群
掃描下方二維碼添加小編微信，小編帶你入伙啦，大牛如雲，讓交流變得簡單。（記得備注姓名-研究方向-單位，防止不小心忽略）