高手题目一答案：nature高级组合图形绘制教程
2020-04-20 本文已影响0人
wentaomicro
写在前面
开学了，大量工作来袭。我在成长，希望可以帮助大家一起成长。按照咱们的约定，在看数目超过50，我应该放出高手题目一的代码。
如果需要我重复高手题目二，请继续点赞，达到100个在看，献上高手题目二代码。点击这里查看
image
准备包和数据
使用内置数据，所以别再问我要了。各位大哥大姐。
library(ggplot2)# 用于绘图library(psych)# 用于相关计算dim(mtcars)
准备参数
由于矩形内是相关，所以需要指定两个阈值。由于是矩形展示，所以需要提供矩形长宽，这里我为了简化，直接作为正方形，所以只需要指定边长。
r.threshold=0.9p.threshold=0.05# 设置矩形半径r = 5#--坐标轴标度长度leng = 0.2
相关计算
#转置mtcars2 = t(mtcars)occor = corr.test(mtcars2,use="pairwise",method="pearson",adjust="fdr",alpha=.05)# occor = cor(t(network_sub),use="pairwise",method="spearman")print("over")occor.r = occor$r # 取相关性矩阵R值occor.p = occor$p # 取相关性矩阵p值occor.r[occor.p > p.threshold|abs(occor.r)<r.threshold] = 0head(occor.r)
矩形展示需要四个分组
模拟四个分组
#---模拟四个分组netClu = data.frame(ID = colnames(mtcars2),group =rep(1:4,length(colnames(mtcars2)))[1:length(colnames(mtcars2))] )netClu$group = as.factor(netClu$group)head(netClu)row.names(netClu) = netClu$ID
> head(netClu) ID group1 Mazda RX4 12 Mazda RX4 Wag 23 Datsun 710 34 Hornet 4 Drive 45 Hornet Sportabout 16 Valiant 2
添加展示丰度柱状图值
netClu = merge(netClu,as.data.frame(colSums(mtcars2)),by = "row.names",all= F)colnames(netClu)[4] = "abu"netClu$abu = netClu$abu/sum(netClu$abu)*30
head(netClu) Row.names ID group abu1 AMC Javelin AMC Javelin 3 1.08896362 Cadillac Fleetwood Cadillac Fleetwood 3 1.56767203 Camaro Z28 Camaro Z28 4 1.39062684 Chrysler Imperial Chrysler Imperial 1 1.56150735 Datsun 710 Datsun 710 3 0.55854886 Dodge Challenger Dodge Challenger 2 1.1181519
四条边分别准备数据
##提取一条边数据gr = dplyr::filter(netClu, group == as.vector(as.factor(names(table(netClu$group)[1]))[1]) )# 根据该组内元素多少准备坐标x1 = seq(from=0, to=r, length.out=table(netClu$group)[1])#--相关点坐标数据构建axis1 = data.frame(row.names = gr$ID,x = x1,y = 0,element = gr$ID)#--丰度坐标构建abu1 = data.frame(row.names = gr$ID,x = x1,y = -gr$abu,element = gr$ID )
其他三个边坐标构造
gr = dplyr::filter(netClu, group == as.vector(as.factor(names(table(netClu$group)[3]))[1]) )x2 = seq(from=0, to=r, length.out=table(netClu$group)[3])axis2 = data.frame(row.names = gr$ID,x = x2,y = r +2,element = gr$ID)abu2 = data.frame(row.names = gr$ID,x = x2,y = gr$abu + r +2 ,element = gr$ID )gr = dplyr::filter(netClu, group == as.vector(as.factor(names(table(netClu$group)[2]))[1]) )y1 = seq(from=0, to=r, length.out=table(netClu$group)[2])axis3 = data.frame(row.names = gr$ID,x = -1,y = y1 + 1,element = gr$ID)abu3 = data.frame(row.names = gr$ID,x = -1 -gr$abu,y = y1 + 1 ,element = gr$ID )gr = dplyr::filter(netClu, group == as.vector(as.factor(names(table(netClu$group)[4]))[1]) )y2 = seq(from=0, to=r, length.out=table(netClu$group)[4])axis4 = data.frame(row.names = gr$ID,x = r +1,y = y2 +1,element = gr$ID)abu4 = data.frame(row.names = gr$ID,x = r +1 + gr$abu,y = y2 + 1 ,element = gr$ID )
图形坐标标度构建
a = data.frame(x = 6,y = seq(from=0, to=-max(abs(abu1$y)), length.out=3))b = data.frame(x = 6 +leng,y = seq(from=0, to=-max(abs(abu1$y)), length.out=3))a1 = data.frame(x = -1,y = seq(from=r +2, to=max(abs(abu2$y)), length.out=3))a1b1 = data.frame(x = -1- leng,y =seq(from=r +2, to=max(abs(abu2$y)), length.out=3))b1a3 = data.frame(x = seq(from=-1, to=-max(abs(abu3$x)), length.out=3),y = 0)a3b3 = data.frame(x = seq(from=-1, to=-max(abs(abu3$x)), length.out=3),y = 0-leng)b3a4 = data.frame(x = seq(from=6, to=max(abs(abu4$x)), length.out=3),y = r+2)a4b4 = data.frame(x = seq(from=6, to=max(abs(abu4$x)), length.out=3),y = r+2+leng)b4
组合矩形坐标和柱状图坐标
point = rbind(axis1,axis2,axis3,axis4)head(point)adun = rbind(abu1,abu2,abu3,abu4)head(adun)colnames(adun) = c("x2","y2","element2")abunpro = merge(point,adun,by = "row.names",all = F)head(abunpro)
可视化效果
ggplot() +geom_point(data = adun,aes(x =x2,y = y2 )) +geom_point(data = point,aes(x =x,y = y ))ggplot() + geom_segment(aes(x = x, y = y, xend = x2, yend = y2), data = abunpro, size = 5)
image
矩形为网络结构 构造网络边和节点
#--节点坐标匹配cor = occor.rplotcord = point#----使用相关矩阵构建网络--network包构建网络#-----g <- network::network(cor, directed=FALSE)#---转化为0-1的相关矩阵# origm <- network::as.matrix.network.adjacency(g) # get sociomatrixedglist <- network::as.matrix.network.edgelist(g)edglist = as.data.frame(edglist)head(edglist)#t添加权重#---network::set.edge.value(g,"weigt",cor)# g# ?get.edge.attribute# row.names(plotcord) = NULLedglist$weight = as.numeric(network::get.edge.attribute(g,"weigt"))head(edglist)edges <- data.frame(plotcord[edglist[, 1], ], plotcord[edglist[, 2], ])head(edges)edges$weight = as.numeric(network::get.edge.attribute(g,"weigt"))##这里将边权重根据正负分为两类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 = as.factor(aaa)colnames(edges) <- c("X1", "Y1","OTU_1", "X2", "Y2","OTU_2","weight","wei_label")edges$midX <- (edges$X1 + edges$X2)/2edges$midY <- (edges$Y1 + edges$Y2)/2head(edges)#--构造边head(netClu)#--返回第一个在第二个中的位置point[,"group"] = netClu[match(point$element,netClu$ID),][3]nodes = pointhead(nodes)head(abunpro)head(netClu)library(dplyr)abunpro = abunpro %>% inner_join(netClu, by = "Row.names")
角度调整
angle = data.frame(group = as.data.frame(table(netClu$group))$Var,c(90,0,90,0))abunpro = abunpro %>% left_join(angle, by = "group")colnames(abunpro)[dim(abunpro)[2]] = "angle"abunpro$hjust = NULLhjust = data.frame(group = as.data.frame(table(netClu$group))$Var,c(1,1,0,0))hjustabunpro = abunpro %>% left_join(hjust, by = "group")colnames(abunpro)[dim(abunpro)[2]] = "hjust"
出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(wei_label)),color = "green", data = edges, size = 0.5) + geom_segment(aes(x = x, y = y, xend = x2, yend = y2,color = group), data = abunpro, size = 10) + # geom_point(aes(x = x, y = y,fill = group,group = group),pch = 21, data = nodes,size = 3) + geom_text(aes(x= x2, y = y2,label = element2,angle = angle,hjust =hjust), data = abunpro)+ scale_colour_brewer(palette = "Set1") + scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) + # labs( title = paste(layout,"network",sep = "_"))+ # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+ # 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())
image.gif
添加坐标
pnet = pnet + geom_segment(aes(x = c(6), y = c(0), xend = c(6), yend = c(-max(abs(abu1$y)))), size = 1) + geom_segment(aes(x = a$x, y = a$y, xend = b$x, yend = b$y),size = 1) + geom_text(aes(x = b$x, y = b$y, label = abs(round(a$y,2))))+ geom_segment(aes(x = c(-1), y = c(r+2), xend = c(-1), yend = c(max(abs(abu2$y)))), size = 1) + geom_segment(aes(x = a1$x, y = a1$y, xend = b1$x, yend = b1$y),size = 1) + geom_text(aes(x = b1$x, y = b1$y, label = abs(round(a1$y,2))-r-2)) + geom_segment(aes(x = c(-1), y = c(0), xend = c(-max(abs(abu3$x))), yend = c(0)), size = 1) + geom_segment(aes(x = a3$x, y = a3$y, xend = b3$x, yend = b3$y),size = 1) + geom_text(aes(x = b3$x, y = b3$y, label = abs(round(a3$x,2))-1)) + geom_segment(aes(x = c(6), y = c(r+2), xend = c(max(abs(abu4$x))), yend = c(r+2)), size = 1) + geom_segment(aes(x = a4$x, y = a4$y, xend = b4$x, yend = b4$y),size = 1) + geom_text(aes(x = b4$x, y = b4$y, label = abs(round(a4$x,2))-r-1))ggsave("./cs3.pdf",pnet,width =15,height = 10)
image
右侧相关
#-------构造右侧右侧数据和另外一组数据的相关。data2 = data.frame(ID = paste("A",1:8,sep = ""), wei = c(1,-1,1,-1,1,-1,1,-1), x1= axis4$x + 5, y1 = axis4$y[8:1] )head(data2 )data3 = cbind(axis4,data2)head(data3)data3$wei = runif(length(data3$x), 0, 0.5)pnet = pnet + geom_point(aes(x = x + 3,y = y),data = axis4,color = "blue",size = 4) + geom_point(aes(x = x1 ,y = y1),data = data2,color = "blue",size = 4) + geom_segment(aes(x = x +3, y = y, xend = x1, yend = y1,size = wei),data = data3,color = "yellow") + geom_text(aes(x= x1,y = y1,label = ID),data = data2,hjust = -1)ggsave("./cs4.pdf",pnet,width =20,height = 10)pnet
image
欢迎加入微生信生物
image.gif
快来微生信生物
轻松一刻**** ◆ ◆
二师兄的日常
二师兄，何许人！小弟亲师兄也，硕士毕业于2018年，你就看着他，就有数不清的意思。在枯燥乏味的科研生活中有着独特的光芒，让我膜拜。如果你感到无力，请关注二师兄，看看他能带给你多少意思。
image
