# 水稻微生物组时间序列分析精讲1-模式图与主坐标轴分析

## 写在前面

《水稻微生物组时间序列分析》的文章，大家对其中图绘制过程比较感兴趣。一上午收到了超30条留言，累计收到41个小伙伴的留言求精讲。

## 图1. 水稻根系微生物随时间变化吗？

A. 水稻全生育期根系微生物组实验设计模式图(以水稻日本晴和IR24为材料,并分别种植于昌平和上庄两地，CP代表北京昌平农场，SZ代表北京海淀上庄)

B-C. 主坐标轴分析(PCoA)展示水稻微生物组随时间变化，其中微生物群落结构主要在第1/2轴上随时间变化(B)，而不同土壤类型主要在第3轴上明显分开(C)

## 图1A. 模式图

https://v.qq.com/x/page/i0633druhc2.html

## 图1B/C. 主坐标轴分析

### 准备工作

``````# Fig.1 PCoA

# Set working enviroment in Rstudio, select Session - Set working directory - To source file location, default is runing directory
``````

``````# clean enviroment object
rm(list=ls())
``````

``````# load related packages
library("ggplot2")
library("vegan")
``````

``````# Set ggplot2 drawing parameter, such as axis line and text size, lengend and title size, and so on.
main_theme = theme(panel.background=element_blank(),
panel.grid=element_blank(),
axis.line.x=element_line(size=.5, colour="black"),
axis.line.y=element_line(size=.5, colour="black"),
axis.ticks=element_line(color="black"),
axis.text=element_text(color="black", size=7),
legend.position="right",
legend.background=element_blank(),
legend.key=element_blank(),
legend.text= element_text(size=7),
text=element_text(family="sans", size=7))
``````

``````# Public file 1. "design.txt"  Design of experiment

# setting subset design
if (TRUE){
sub_design = subset(design,groupID %in% c("A50Cp0","A50Cp1","A50Cp10","A50Cp112","A50Cp119","A50Cp14","A50Cp2","A50Cp21","A50Cp28","A50Cp3","A50Cp35","A50Cp42","A50Cp49","A50Cp63","A50Cp7","A50Cp70","A50Cp77","A50Cp84","A50Cp91","A50Cp98","A50Sz0","A50Sz1","A50Sz10","A50Sz118","A50Sz13","A50Sz2","A50Sz20","A50Sz27","A50Sz3","A50Sz34","A50Sz41","A50Sz48","A50Sz5","A50Sz56","A50Sz62","A50Sz69","A50Sz7","A50Sz76","A50Sz83","A50Sz90","A50Sz97","HNCp112","HNCp119","HNSz118","IR24Cp0","IR24Cp1","IR24Cp10","IR24Cp112","IR24Cp119","IR24Cp14","IR24Cp2","IR24Cp21","IR24Cp28","IR24Cp3","IR24Cp35","IR24Cp42","IR24Cp49","IR24Cp63","IR24Cp7","IR24Cp70","IR24Cp77","IR24Cp84","IR24Cp91","IR24Cp98","IR24Sz0","IR24Sz1","IR24Sz10","IR24Sz118","IR24Sz13","IR24Sz2","IR24Sz20","IR24Sz27","IR24Sz3","IR24Sz34","IR24Sz41","IR24Sz48","IR24Sz5","IR24Sz56","IR24Sz62","IR24Sz69","IR24Sz7","IR24Sz76","IR24Sz83","IR24Sz90","IR24Sz97","soilCp0","soilCp42","soilCp49","soilCp57","soilCp63","soilCp77","soilCp84","soilCp91","soilCp98","soilSz0","soilSz41","soilSz48","soilSz54","soilSz62","soilSz76","soilSz83","soilSz90","soilSz97") ) # select group1
}else{
sub_design = design
}

print(paste("Number of group: ",length(unique(sub_design\$group)),sep="")) # show group numbers
``````

``````#  PCoA bray_curtis
``````

``````# subset matrix and design
idx = rownames(sub_design) %in% colnames(bray_curtis)
sub_design = sub_design[idx,]
bray_curtis = bray_curtis[rownames(sub_design), rownames(sub_design)] # subset and reorder distance matrix
``````

``````# cmdscale {stats}, Classical multidimensional scaling (MDS) of a data matrix. Also known as principal coordinates analysis
pcoa = cmdscale(bray_curtis, k=4, eig=T) # k is dimension, 3 is recommended; eig is eigenvalues
points = as.data.frame(pcoa\$points) # get coordinate string, format to dataframme
colnames(points) = c("x", "y", "z","a")
eig = pcoa\$eig
``````

``````points = cbind(points, sub_design[match(rownames(points), rownames(sub_design)), ])
``````

``````# plot PCo 1 and 2
p = ggplot(points, aes(x=x, y=y, color=day, shape = compartment))
p = p + geom_point(alpha=.7, size=2) +
labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),
y=paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4), "%)", sep=""),
title="bray_curtis PCoA") + main_theme
``````

``````p
ggsave("beta_pcoa_day_bray_curtis_default.pdf", q, width = 4, height = 2.5)
``````

``````# scale_color_gradientn 按多种颜色连续着色，如彩虹色
# topo.colors(), rainbow(), heat.colors(), terrain.colors(), cm.colors(), RColorBrewer::brewer.pal()
q
ggsave("beta_pcoa_day_bray_curtis_rainbow.pdf", q, width = 4, height = 2.5)
``````

PCoA有众多维度，通常看1/2轴可解析的差异也最大。但有时你研究的目标末必是最大差异，可以进一步往下找，如1/3，1/4，3/4轴。

``````# plot PCo 1 and 3
points\$siteXcompt=paste(points\$site,points\$compartment,sep = "")

p = ggplot(points, aes(x=x, y=z, color=day, shape = siteXcompt))
p = p + geom_point(alpha=.7, size=2) +
labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),
y=paste("PCoA 3 (", format(100 * eig[3] / sum(eig), digits=4), "%)", sep=""),
title="bray_curtis PCoA") + main_theme
p
ggsave("beta_pcoa_day_bray_curtis3.pdf", p, width = 4, height = 2.5)
``````

Citition:
Zhang, J., Zhang, N., Liu, Y.X., Zhang, X., Hu, B., Qin, Y., Xu, H., Wang, H., Guo, X., Qian, J., et al. (2018). Root microbiota shift in rice correlates
with resident time in the field and developmental stage. Sci China Life Sci 61, https://doi.org/10.1007/s11427-018-9284-4

## 写在后面

https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA