写在前面

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

我们也争取花时间把此文的原始代码整理并精讲，祝有需要的小伙伴能有收获。

本系列按原文4幅组图，共分为4节。本文是第一节，模式图与主坐标轴分析。

先回顾一下图1的内容。

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

实验设计与群落整体结构(主坐标轴PCoA)分析

图1. 田间水稻微生物组随生育时间变化

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

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

方法说明：图1A采用Aodbe Illustrator手绘，具体方法可参考youbute上使用AI绘制花的教程，https://www.youtube.com/watch?v=WSxemBP-gZQ。
图1B/C是基于Bray-Curtis距离进行的PCoA分析，采用散点图展示，并按时间顺序填充彩虹色，按不同compartment设置形状。图1B展示PCo1/2轴，组间最大差异为不同compartment与时间梯度上的变化。图1C展示PCo1/3轴，可进一步看到1轴的差异与时间变化一致，而3轴可以很好分开不同地点。

图1A. 模式图

绘制模式图，工具有非常多的选择。有的大神用PPT可以画出美到极致的模式图，有的人用喜欢用PS。在这里我们推荐用AI，因为它是学术界矢量图绘制排版神器，没有之一。

AI是Adobe Illustrator软件的缩写，对于各种软件生成的矢量图，混合编辑，轻松满足各类杂志的所有要求。没有基础的小伙伴可以学习下文：

• 教师节献礼 - 文章用图的修改和排版

本文中是用AI绘制了水稻不同生育期的矢量图，绘制还是要求有一定基础和思路，以及AI的基本使用技巧才能完成。

这里推荐来自youbute上使用AI绘制花和蝴蝶的教程，https://www.youtube.com/watch?v=WSxemBP-gZQ。

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

个人感觉软件操作比较容易学习，难点是如何设计合理的模式图，及基本的绘图技术。设计、构思、修改都是极花时间的，只有原作者对自己的项目理解最深刻，最容易构思较好的模式图，具体的绘制建议拿草图找会画画的同学合作。

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

主坐标轴分析需要两个输入文件，一个是实验设计，即样本与分组(时间)对应表。另一个是样本间的距离矩阵(距离矩阵可通过OTUs直接计算获得，详见《扩增子分析流程》中第6节)。

下面开始数据筛选与绘图，请按文末说明下载数据和代码，在Rstudio中即可重现。

准备工作

打开代码文件，设置工作目录

# Fig.1 PCoA

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

每次绘图前(list=ls())清空工作环境，防止有之前的旧变量导致绘图结果中有错误

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

加载本分析需要的包，主要包括生态统计的vegan包和一年引用超7000次的绘图神器ggplot2，没有安装过此包的朋友请使用Rstudio中Packages页面手动安装

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

设置图形输出的基本样式，即theme，每个人都有自己的风格，这是我的习惯风格，如字体7号是Nature杂志推荐的字号。经验积累一个自己的theme，可以为后期AI排版节约很多时间，也有自己的风格。

# 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))

读取实验设计并筛选

现在的实验一般都有特别多的组，如本研究包括两品种、两地点，再乘以时间点近100组，500个样品。肯定会有人为造成的异常样本，如何通过实验分析结果、配合实验记录排除人为错误或不可控因素造成的影响呢？那就是仔细观察数据和实验记录，综合判断后筛选样本，是十分必要的。

我们的分析中也发现了两个时间点样本的异常，根据农场的记录对应，发现了这两个点分别进行了大规模的施肥和除虫。将确实对分析影响较大的异常样本和组剔除后再次分析。

# Public file 1. "design.txt"  Design of experiment
design = read.table("../data/design.txt", header=T, row.names= 1, sep="\t") 

# 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
bray_curtis = read.table("../data/bray_curtis_otu_table_css.txt", sep="\t", header=T, check.names=F)

距离矩阵与实验设计的交叉筛选

这步是非常必要的，因为实验设计在不断的分析中，会删除一些异常样本。而OTU表中也会有些样本数据量过低、过高而被删除。大量样本时不完全对应，需要进行交叉筛选保持实验设计与数据矩阵一致。

# 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

结果有很多信息，主要提取结果中的坐标位置和各轴的解析差异的比例

# 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

添加样品组信息：合并PCoA坐标与实验设计

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

绘图：散点图展示样品坐标，XY轴标签显示可解析差异的比例，按时间序列为分组连序着色，按不同取样部分显示不同形状

# 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

预览结果并保存：注意图片输出PDF格式，可以在AI中进一步编辑文件和线条，图片大小4 x 2.5适合大多数杂志的1/2栏大小

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

这个图是把时间变化展示出来了，但默认的深蓝到浅蓝，太低调，不够妖艳。我个人喜欢R ggplot2绘图另一个原因是喜欢它的彩虹色系统。

彩虹色绘制第1/2主轴的时间变化规律

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

现在看时间变化是不是清楚多了，纯个人感觉，也行有的朋友不喜欢这种风格，更有一些杂志社严格要求不允许使用彩虹色。那就要看你根据杂志要求和自己的表现目标调整了，上面代码注释行中也提了众多可调整的函数方案可选。

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

# 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) +
  scale_color_gradientn(colours=rainbow(7)) +
  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)

采用1/3轴组合，即展示出时间序列变化，又展示出了两地点昌平(Cp)和上庄(Sz)间的明显差异。

本分析的全部文件和代码，会在 https://github.com/YongxinLiu/RiceTimeCourse 上持续更新，也可以后台回复“时间序列”获得百度云下载链接。

如果本文分享的技术帮助了你的科研，欢迎引用下文，支持国产SCI越来越好。

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

点我下载PDF

猜你喜欢

• 10000+：肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊
• 系列教程：微生物组入门 Biostar 微生物组  宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索  Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记  云协作 公众号
• 编程模板 Shell  R Perl
• 生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外1400+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA