查看原文
其他

扩增子统计绘图2散点图:Beta多样性

2017-08-22 刘永鑫 宏基因组

点击上方蓝色「宏基因组」关注我们!专业干货每日推送!

写在前面

优秀的作品都有三部分曲,如骇客帝国、教父、指环王等。

扩增子系列课程也分为三部曲:

第一部《扩增子图表解读》:加速大家对同行文章的解读能力。

第二部《扩增子分析解读》:学习数据分析的基本思路和流程。

第三部《扩增子统计绘图》:即是对结果进行可视和统计检验,达到出版级的图表结果。

《扩增子统计绘图》系列文章介绍

《扩增子统计绘图》是之前发布的《扩增子图表解读》和《扩增子分析解读》的进阶篇,是在大家可以看懂文献图表,并能开展标准扩增子分析的基础上,进行结果的统计与可视化。其章节设计与《扩增子图表解读》对应,为八节课八种常用图形(箱线图、散点图、热图、曼哈顿图、火山图、维恩图、三元图和网络图),基本满足文章常用的图片种类需求。

也适合对公司标准化分析返回结果的进一步统计、可视化及美化,达到出版级别,冲击高分文章。

本部分练习所需文件位于百度网盘,链接:http://pan.baidu.com/s/1hs1PXcw 密码:y33d。

绘制Beta多样性线散点图

绘图和统计全部为R语言,建议复制代码,在Rstuido中运行,并设置工作目录为存储之前分析结果文件的result目录。

主坐标轴分析(PCoA)展示所有样品间的最大差异

采用PCoA展示样品间的最大差异,代码以bray_curtis为例,其它距离只需替换为weighted_unifrac,unweighted_unifrac即可。

# 运行前,请在Rstudio中菜单栏选择“Session - Set work directory -- Choose directory”,弹窗选择之前分析目录中的result文件夹 # 安装相关软件包,如果末安装改为TRUE运行即可安装 if (FALSE){  source("https://bioconductor.org/biocLite.R")  biocLite(c("ggplot2","vegan")) } # 加载相关软件包 library("ggplot2") library("vegan") # 读入实验设计 design = read.table("design.txt", header=T, row.names= 1, sep="\t") # 读入bray_curtis的距离矩阵 bray_curtis = read.table("beta/bray_curtis_otu_table_css.txt", sep="\t", header=T, check.names=F) # 过滤数据并排序 idx = rownames(design) %in% colnames(bray_curtis) sub_design = design[idx,] bray_curtis = bray_curtis[rownames(sub_design), rownames(sub_design)] # subset and reorder distance matrix # 将距离矩阵进行主坐标轴分析 pcoa = cmdscale(bray_curtis, k=3, 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") eig = pcoa$eig points = cbind(points, sub_design[match(rownames(points), rownames(sub_design)), ]) # 绘制主标准轴的第1,2轴 p = ggplot(points, aes(x=x, y=y, color=genotype)) +  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") p ggsave("beta_pcoa_bray_curtis.pdf", p, width = 5, height = 3) ggsave("beta_pcoa_bray_curtis.png", p, width = 5, height = 3)

详细的图片讲解,可参考2散点图:Beta多样性,组间整体差异分析
图中WT和OE在第一轴明显分开的,但WT与KO间区分不明显,是否存在显著差别呢。

我们以WT和OE为例,统计组间是否存在显著差异。

# 统计WT与KO间是否有差异显著 # 取两组实验设计子集 design2 = subset(sub_design, genotype %in% c("WT","KO")) # 获取对应的子距离矩阵并排序 sub_dis_table = bray_curtis[rownames(design2),rownames(design2)] # 计算距离矩阵 sub_dis_table <- as.dist(sub_dis_table, diag = FALSE, upper = FALSE) # 统计按genotype分组下,组间差异的显著性水平;检验10000次 adonis_table = adonis(sub_dis_table~genotype, data=design2, permutations = 10000) # 获得pvalue值 adonis_pvalue = adonis_table$aov.tab$`Pr(>F)`[1] # 显示组间的pvalue值 adonis_pvalue

pvalue值结果不是一个确定的值,是多次统计的结果,但差别不大,每次大约为0.01左右,即WT与OE存在显著差异,但还是不极显著。读者可以自行计算WT与OE试试,看看P值有多显著。

限制条件主坐标轴分析(CCA)展示组间的最大差异

通常PCoA展示的是所有样品间的最大差异,而实验中想要找的是组间的差异,就需要限制条件的主坐标轴分析。

# CCA分析功能函数 variability_table = function(cca){  chi = c(cca$tot.chi, cca$CCA$tot.chi, cca$CA$tot.chi)  variability_table = cbind(chi, chi/chi[1])  colnames(variability_table) = c("inertia", "proportion")  rownames(variability_table) = c("total", "constrained", "unconstrained")  return(variability_table) } # 读入CSS标准化的OTU表,并与实验设计比对筛选和数据重排 otu_table = read.table("otu_table_css.txt", sep="\t", header=T, row.names= 1) # CSS norm otu table idx = rownames(sub_design) %in% colnames(otu_table) sub_design = sub_design[idx,] sub_otu_table = otu_table[, rownames(sub_design)] # Constrained analysis OTU table by genotype capscale.gen = capscale(t(sub_otu_table) ~ genotype, data=sub_design, add=F, sqrt.dist=T, distance="bray") # ANOVA-like permutation analysis perm_anova.gen = anova.cca(capscale.gen) # generate variability tables and calculate confidence intervals for the variance var_tbl.gen = variability_table(capscale.gen) eig = capscale.gen$CCA$eig variance = var_tbl.gen["constrained", "proportion"] p.val = perm_anova.gen[1, 4] # extract the weighted average (sample) scores points = capscale.gen$CCA$wa[, 1:2] points = as.data.frame(points) colnames(points) = c("x", "y") points = cbind(points, sub_design[match(rownames(points), rownames(sub_design)),]) # plot CPCo 1 and 2 p = ggplot(points, aes(x=x, y=y, color=genotype)) +  geom_point(alpha=.7, size=1.5) +  labs(x=paste("CPCoA 1 (", format(100 * eig[1] / sum(eig), digits=4), "%)", sep=""),       y=paste("CPCoA 2 (", format(100 * eig[2] / sum(eig), digits=4), "%)", sep="")) +  ggtitle(paste(format(100 * variance, digits=3), " % of variance; p=",format(p.val, digits=2),sep="")) p ggsave(paste( "CPCoA.pdf", sep=""), p, width = 5, height = 3) ggsave(paste( "CPCoA.png", sep=""), p, width = 5, height = 3)

图中三个组能明显分开,代表组间存在一致的差异。顶部展示21.2%表示组间的差异占总体的比例,p=0.001表示组间有显著差异。两轴百分比是此平面下可解释差异的百分比。
详细的图片讲解,可参考2散点图:Beta多样性,组间整体差异分析  

想了解更多宏基因组、16S分析相关文章,

快关注“宏基因组”公众号,干货第一时间推送。

系统学习生物信息,快关注“生信宝典”,

那里有几千志同道合的小伙伴一起学习。

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

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