查看原文
其他

水稻微生物组时间序列分析2a-相关分析

LiuYX 宏基因组 2019-07-10

写在前面

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

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

本系列按原文4幅组图,共分为4节。本文是第二节a,相关分析。图2的散点图和拟合部分内容独立,将在下节2b分享。

前情提要

先了解一下图2的内容。

图2. 微生物组随时间变化的规律

图2. 田间水稻根系微生物组在8-10周趋于稳定。

A-D. 对两个水稻品种分别在两地进行的连续微生物组调查结果相关分析,发现8-10周后群落结构趋于稳定。

E. 所有时间点距离水稻最后取样点的Bray-Curtis距离。发现土壤呈现小幅波动变化,而根系呈现出先快后慢,逐渐趋近的变化规律。

F. 不同水稻品种在两个地点间的距离变化,发现土壤差异稳定,而根系微生物组差异随时间增长而趋于一致。

方法简介:A-D采用R的cor()函数计算pearson相关系数,并使用Corrplot包展示,时间轴使用pheatmap绘制热图展示。

E-F基于vegan包计算的所有样品两两间Bray-Curtis距离。分别挑选距离终点的距离,和两地间的距离与时间序列上的关系,并采用ggplot2可视化散点图,并添加拟合曲线方便观察变化规律。

图2A-D. 相关性corrplot

A-D为四个相当类型图,只是分别两个地点的两个品种进行分析,进一步説在不同地点和不同品种的变量下仍然存在稳定的变化规律。这里仅以图2A在北京种植的日本晴水稻品种为例进行代码说明。

输入文件只有实验设计和OTU表,具体文件见文末说明。

清空工作环境和加载包

# Set working enviroment in Rstudio, select Session - Set working directory - To source file location, default is runing directory rm(list=ls()) # clean enviroment object library("corrplot") library("pheatmap") library(ggcorrplot)

读入实验设计和OTU表

# Public file 1. "design.txt"  Design of experiment design = read.table("../data/design.txt", header=T, row.names= 1, sep="\t") # Public file 2. "otu_table.txt"  raw reads count of each OTU in each sample otu_table = read.delim("../data/otu_table.txt", row.names= 1,  header=T, sep="\t")

实验设计筛选:这里只筛选昌平的日本晴相关组(A50)

# setting subset design if (TRUE){    sub_design = subset(design,groupID %in% c("A50Cp0","A50Cp1","A50Cp2","A50Cp3","A50Cp7","A50Cp10","A50Cp14","A50Cp21","A50Cp28","A50Cp35","A50Cp42","A50Cp49","A50Cp63","A50Cp70","A50Cp77","A50Cp84","A50Cp91","A50Cp98","A50Cp112","A50Cp119") ) # select group1 }else{    sub_design = design } sub_design$group=sub_design$groupID

设置组顺序,这是必须的,否则时间10,100天会位于1的后面。

# Set group order if ("TRUE" == "TRUE") {    sub_design$group  = factor(sub_design$group, levels=c("A50Cp0","A50Cp1","A50Cp2","A50Cp3","A50Cp7","A50Cp10","A50Cp14","A50Cp21","A50Cp28","A50Cp35","A50Cp42","A50Cp49","A50Cp63","A50Cp70","A50Cp77","A50Cp84","A50Cp91","A50Cp98","A50Cp112","A50Cp119"))    }else{sub_design$group  = as.factor(sub_design$group)} print(paste("Number of group: ",length(unique(sub_design$group)),sep="")) # show group numbers

筛选后的实验设计样本与OTU表交叉筛选

# sub and reorder subdesign and otu_table idx = rownames(sub_design) %in% colnames(otu_table) sub_design = sub_design[idx,] count = otu_table[, rownames(sub_design)]

OTU表标准化为百分比,在R中只需一小行代码

norm = t(t(count)/colSums(count,na=T)) * 100 # normalization to total 100

按组合并:因为样本太多,一小部分过百个样本,展示太乱,按组求均值,组间比较更容易发现规律

# Pearson correlation among groups sampFile = as.data.frame(sub_design$group,row.names = row.names(sub_design)) colnames(sampFile)[1] = "group" mat = norm mat_t = t(mat) mat_t2 = merge(sampFile, mat_t, by="row.names") mat_t2 = mat_t2[,-1] mat_mean = aggregate(mat_t2[,-1], by=mat_t2[1], FUN=mean) # mean mat_mean_final = do.call(rbind, mat_mean)[-1,] geno = mat_mean$group colnames(mat_mean_final) = geno

计算相关系数,并保留3位小数

sim=cor(mat_mean_final,method="pearson") sim=round(sim,3)

先使用ggcorrplot绘制相关矩阵

pdf(file="ggcorrplot_pearson_A50Cp.pdf", height = 2.5, width = 2.5) ggcorrplot(sim, type = "upper", colors = c("green", "yellow", "red")) # , method = "circle" dev.off()


这是ggcorrplot绘制的默认效果,还是很不错的。
详细学习见 R相关矩阵可视化包ggcorrplot

但我个人更喜欢更老的corrplot包。

# 人为设置颜色度 col1 <- colorRampPalette(c("green", "green", "red")) pdf(file="corplot_pie_pearson_A50Cp.pdf", height = 2.5, width = 2.5) corrplot(sim, method="pie", type="lower", col=col1(100)) # , diag=F , na.label = "1" dev.off()


采用corrplot的饼形图样式,颜色为绿至红。

更多教程见:R画月亮阴晴圆缺:corrplot绘图相关系数矩阵

生成图例

col1 <- colorRampPalette(c("green", "red")) corrplot(sim, method="pie", type="lower", col=col1(100)) # , diag=F , na.label = "1" # 生成时间热图,分别为土和植物的 time1 = c(0,1,2,3,7,10,14,21,28,35,42,49,63,70,77,84,91,98,112,119) time2 = c(0,41,48,54,62,77,84,90,97,119,0,0,0,0,0,0,0,0,0,0) time=data.frame(time1,time2) pheatmap(time, cluster_rows = F,  cluster_cols = F) pheatmap(time, cluster_rows = F,  cluster_cols = F, filename = "corplot_pie_legend_time.pdf" ,width=2, height=4)

植物和土壤的时间热图,我们采用pheatmap绘制,再用用AI添加在图中的。

本分析的全部文件和代码,会在 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

猜你喜欢

10000+:肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:生信宝典 学术图表 高分文章 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树

必备技能:提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具:16S预测培养基 生信绘图

科研经验:云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

猜你喜欢

10000+:肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大  Cell微生物专刊

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:生信宝典 学术图表 高分文章 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树

必备技能:提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具:16S预测培养基 生信绘图

科研经验:云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

猜你喜欢

10000+:肠道细菌 人体上的生命 宝宝与猫狗 梅毒狂想曲 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊

系列教程:微生物组入门 Biostar 微生物组  宏基因组

专业技能:生信宝典 学术图表 高分文章 不可或缺的人

一文读懂:宏基因组 寄生虫益处 进化树

必备技能:提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析:图表解读 分析流程 统计绘图

16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

在线工具:16S预测培养基 生信绘图

科研经验:云笔记  云协作 公众号

编程模板 Shell  R Perl

生物科普  生命大跃进  细胞暗战 人体奥秘  

写在后面

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

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

点击阅读原文,跳转最新文章目录阅读

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存