扩增子统计绘图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 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文献阅读和分析相关文章，快关注“宏基因组”公众号，干货第一时间推送。

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