本文从前文 MicEco:计算Sloan随机性的另一方法 得到的结果中继续画图,希望得到最开始mSystems中类似的图。
1#随便找了一个OTU表
2otu = read.table(file = "otu.txt",sep="\t",row.names=1,header=T)
3
4library(MicEco)
5library(ggplot2)
6res = neutral.fit(t(otu))
7
8##########画图
9m = res[[1]][1]
10r2 = res[[1]][3]
11out = res[[2]]
先画三条线,预测的发生率及上下界
1p1 = ggplot() +
2 geom_line(data = out,aes(x=log(p),y=freq.pred),size = 1.2,linetype = 1)+
3 geom_line(data = out,aes(x=log(p),y=Lower),size = 1.2,linetype = 2)+
4 geom_line(data = out,aes(x=log(p),y=Upper),size = 1.2,linetype = 2)+
5 xlab("log10(mean relative abundance)")+ylab("Occurrence frequency")
6p1
加散点图
1#判断点的位置,分成三组,红-绿-黄
2require(plyr)
3out2 = mutate(out, group=NA)
4out2$group[out[,2]<out[,4]]="yellow" ##低于下界
5out2$group[out[,2]>out[,5]]="red" ##高于上界
6out2$group[(out[,2]>=out[,4])&(out[,2]<=out[,5])]="green"##中间
7
8mycols<-c("green","red1","orange")
9p2 = p1 + geom_point(data = out2,aes(x=log(p),y=freq,color = group),size = 2)+
10 scale_colour_manual(values = mycols)+
11 annotate("text",x=-12.5,y=0.25,label=paste("m = ",round(m,3),sep=''),size=7)+
12 annotate("text",x=-12.5,y=0.3,label=paste("R2 = ",round(r2,3),sep=''),size=7)
13p2
格式调整
1plot_theme = theme(panel.background=element_blank(),
2 panel.grid=element_blank(),
3 axis.line.x=element_line(size=.5, colour="black"),
4 axis.line.y=element_line(size=.5, colour="black"),
5 axis.ticks=element_line(color="black"),
6 axis.text=element_text(color="black", size=24),
7 legend.position="none", ##right
8 legend.background=element_blank(),
9 legend.key=element_blank(),
10 legend.text= element_text(size=24),
11 text=element_text(family="sans", size=24)
12)
13p3 = p2 + plot_theme;p3
计算三组OTU的比例并画饼图
1low = nrow(out2[out2[,6]== "yellow",])
2med = nrow(out2[out2[,6]== "green",])
3high = nrow(out2[out2[,6]== "red",])
4
5type <- c('med','high','low')
6nums <- c(med,high,low)
7df <- data.frame(type = type, nums = nums)
8label_value <- paste('', round(df$nums/sum(df$nums) * 100, 1), '%', sep = '')
9label_value
10label <- paste(df$type, label_value, sep = ' ')
11label
12p4 <- ggplot(data = df, aes(x = 1, y = nums, fill = type)) +
13 geom_bar(stat = 'identity', position = 'stack', width = 0.5)+
14 scale_fill_manual(name='',
15 labels=c(label[2], label[3], label[1]),
16 values=c("red1","orange","green"))+
17 coord_polar(theta = 'y')+
18 labs(x = '', y = '', title = '')+
19 theme(axis.text = element_blank())+
20 theme(axis.ticks = element_blank())+
21 theme(legend.position = "right")
22
24p4
25p4 = p4+theme(panel.background=element_blank(),
26 panel.grid=element_blank(),
27 legend.background=element_blank(),
28 legend.key=element_blank(),
29 legend.text= element_text(size=24))
30
31p4
拼图
1require(ggimage)
2g <- p3 + geom_subview(subview = p4 + theme_void(), x=-11.5, y=0.75, w=5, h=5)
3g
4#dev.off()
现在已经比较像了。当然还可以继续美化~