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

版权声明:本文为“宏基因组”公众号原创,未经博主允许不得转载。 https://blog.csdn.net/woodcorpse/article/details/77458783

写在前面

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

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

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

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

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

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

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

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

本部分练习所需文件位于百度网盘,链接: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)

image
详细的图片讲解,可参考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)

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

想了解更多宏基因组、16S文献阅读和分析相关文章,快关注“宏基因组”公众号,干货第一时间推送。
image

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

展开阅读全文

没有更多推荐了,返回首页