写在前面我们因此计划花时间把此文的原始代码整理并精讲，祝有需要的小伙伴能有收获。 本系列按原文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./watch?v=WSxemBP-gZQ。 图1A. 模式图绘制模式图，工具有非常多的选择。有的大神用PPT可以画出美到极致的模式图，有的人用喜欢用PS。在这里我们推荐用AI，因为它是学术界矢量图绘制排版神器，没有之一。 AI是Adobe Illustrator软件的缩写，对于各种软件生成的矢量图，混合编辑，轻松满足各类杂志的所有要求。没有基础的小伙伴可以学习下文： 本文中是用AI绘制了水稻不同生育期的矢量图，绘制还是要求有一定基础和思路，以及AI的基本使用技巧才能完成。 这里推荐来自youbute上使用AI绘制花和蝴蝶的教程，https://www./watch?v=WSxemBP-gZQ。
图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 objectrm(list=ls()) 加载本分析需要的包，主要包括生态统计的vegan包和一年引用超7000次的绘图神器ggplot2，没有安装过此包的朋友请使用Rstudio中Packages页面手动安装 # load related packageslibrary('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 experimentdesign = read.table('../data/design.txt', header=T, row.names= 1, sep='\t') # setting subset designif (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_curtisbray_curtis = read.table('../data/bray_curtis_otu_table_css.txt', sep='\t', header=T, check.names=F)距离矩阵与实验设计的交叉筛选 这步是非常必要的，因为实验设计在不断的分析中，会删除一些异常样本。而OTU表中也会有些样本数据量过低、过高而被删除。大量样本时不完全对应，需要进行交叉筛选保持实验设计与数据矩阵一致。 # subset matrix and designidx = 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 analysispcoa = cmdscale(bray_curtis, k=4, eig=T) # k is dimension, 3 is recommended; eig is eigenvaluespoints = as.data.frame(pcoa$points) # get coordinate string, format to dataframmecolnames(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 2p = 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栏大小 pggsave('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))qggsave('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 3points$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_themepggsave('beta_pcoa_day_bray_curtis3.pdf', p, width = 4, height = 2.5)采用1/3轴组合，即展示出时间序列变化，又展示出了两地点昌平(Cp)和上庄(Sz)间的明显差异。 |
|
|
