|
上周五我们分享了3月底发表的的
《水稻微生物组时间序列分析》的文章,大家对其中图绘制过程比较感兴趣。一上午收到了超30条留言,累计收到41个小伙伴的留言求精讲。
我们也争取花时间把此文的原始代码整理并精讲,祝有需要的小伙伴能有收获。
本系列按原文4幅组图,共分为4节。本文是第一节,模式图与主坐标轴分析。
先回顾一下图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轴可以很好分开不同地点。
绘制模式图,工具有非常多的选择。有的大神用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
个人感觉软件操作比较容易学习,难点是如何设计合理的模式图,及基本的绘图技术。设计、构思、修改都是极花时间的,只有原作者对自己的项目理解最深刻,最容易构思较好的模式图,具体的绘制建议拿草图找会画画的同学合作。
主坐标轴分析需要两个输入文件,一个是实验设计,即样本与分组(时间)对应表。另一个是样本间的距离矩阵(距离矩阵可通过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
为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外1400+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
学习扩增子、宏基因组科研思路和分析实战,关注“宏基因组”
点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 16:26
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社