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

写在前面

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

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

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

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

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

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

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

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

本部分练习所需文件位于百度网盘，链接：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 eigenvaluespoints = as.data.frame(pcoa$points) # get coordinate string, format to dataframmecolnames(points) = c("x", "y", "z") eig = pcoa$eigpoints = 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")pggsave("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 tableidx = 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 genotypecapscale.gen = capscale(t(sub_otu_table) ~ genotype, data=sub_design, add=F, sqrt.dist=T, distance="bray") # ANOVA-like permutation analysisperm_anova.gen = anova.cca(capscale.gen)# generate variability tables and calculate confidence intervals for the variancevar_tbl.gen = variability_table(capscale.gen)eig = capscale.gen$CCA$eigvariance = var_tbl.gen["constrained", "proportion"]p.val = perm_anova.gen[1, 4]# extract the weighted average (sample) scorespoints = 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 2p = 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="")) pggsave(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文献阅读和分析相关文章，快关注“宏基因组”公众号，干货第一时间推送。

系统学习生物信息，快关注“生信宝典”，那里有几千志同道合的小伙伴一起学习。