ggplot2版聚类物种丰度堆叠图
文章目录
- 写在前面
- 加载依赖关系
- 导入数据
- ggtree绘制聚类树
- 物种组成数据
- 整理成facet需要的格式
- 保证颜色填充独立性
- 分面组合树和柱图
- 修改配色
- ggtree调整布局
- 添加样本其他信息
- 树+柱+堆叠图组合
- 猜你喜欢
- 写在后面
写在前面
随着研究的逐渐深入,我们对绘图的要求越来越高,各种之前使用的较少的图形如今追求热度和新颖程度,都开始逐渐在大文章中显现。如下图。
这是最近刚发表于Nature Ecology & Evolution中的图1b。如何绘制呢?
Thorsten Thiergart, Paloma Durán, Thomas Ellis, Nathan Vannier, Ruben Garrido-Oter, Eric Kemen, Fabrice Roux, Carlos Alonso-Blanco, Jon Ågren, Paul Schulze-Lefert & Stéphane Hacquard. Root microbiota assembly and adaptive differentiation among European Arabidopsis populations. Nature Ecology & Evolution 4, 122-131, doi:10.1038/s41559-019-1063-3 (2020).
这次的聚类加物种丰度展示让我们学习一波。之前推出了用R语言的plot绘制的教程。
- R语言绘制带聚类树的堆叠柱形图
但修改细节仍比较麻烦。今天更新基于ggplot2系统的教程。
加载依赖关系
这里的ggtree需要使用19年7月以后的版本,因为这以后的版本才支持将聚类结果转化为树结构。
如果你的Bioconductor版本较旧,可能一直会安装旧版ggtree。升新方法如下:
## 先卸载先前的安装控制程序
remove.packages(c("BiocInstaller", "BiocManager", "BiocVersion"))## 再安装新版程序
install.packages("BiocManager")
BiocManager::install(update=TRUE, ask=FALSE)
library("ggplot2")
library("ggdendro")
# library(remotes)
library(phyloseq)
library(tidyverse)
library(ggtree)
library( ggstance)
# library(amplicon)
vegan_otu = function(physeq){OTU = otu_table(physeq)if(taxa_are_rows(OTU)){OTU = t(OTU)}return(as(OTU,"matrix"))
}vegan_tax <- function(physeq){tax <- tax_table(physeq)return(as(tax,"matrix"))
}
导入数据
# 从R数据文件中读入
# ps = readRDS("data/ps_liu.rds")# 从文件读取
metadata = read.table("http://210.75.224.110/github/EasyAmplicon/data/metadata.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F)
otutab = read.table("http://210.75.224.110/github/EasyAmplicon/data/otutab.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F)
taxonomy = read.table("http://210.75.224.110/github/EasyAmplicon/data/taxonomy.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F)# 提取两个表中共有的ID
# Extract only those ID in common between the two tables
idx = rownames(otutab) %in% rownames(taxonomy)
otutab = otutab[idx,]
taxonomy = taxonomy[rownames(otutab),]# 使用amplicon包内置数据
# data("metadata")
# data(otutab)# 导入phyloseq(ps)对象
ps = phyloseq(sample_data(metadata),otu_table(as.matrix(otutab), taxa_are_rows=TRUE), tax_table(as.matrix(taxonomy)))
ggtree绘制聚类树
# 样本间距离类型:Bray-Curtis
dist = "bray"
# phyloseq(ps)对象标准化
ps1_rela = transform_sample_counts(ps, function(x) x / sum(x) )
# 导出OTU表
otu = as.data.frame(t(vegan_otu(ps1_rela)))
# 预览
otu[1:3,1:3]
#计算距离矩阵
unif = phyloseq::distance(ps1_rela , method=dist)
# 聚类树,method默认为complete
hc <- hclust(unif, method = "complete")
# 对树分组
clus <- cutree(hc, 3)
# 提取树中分组的标签和分组编号
d = data.frame(label = names(clus), member = factor(clus))
# 提取样本元数据
map = as.data.frame(sample_data(ps))
# 合并树信息到样本元数据
dd = merge(d,map,by = "row.names",all = F)
row.names(dd) = dd$Row.names
dd$Row.names = NULL
dd[1:3,1:3]# ggtree绘图 #----
p = ggtree(hc) %<+% dd + geom_tippoint(size=5, shape=21, aes(fill=factor(Group), x=x)) + # geom_tiplab(aes(label=Group), size=3, hjust=.5) +geom_tiplab(aes(color = Group,x=x*1.2), hjust=1)# theme_dendrogram(plot.margin=margin(6,6,80,6))# 这是聚类图形的layout
p
物种组成数据
# 指定物种组成的选项
i = ps # 指定输入数据
j = "Phylum" # 使用门水平绘制丰度图表
rep = 6 # 重复数量是6个
Top = 10 # 提取丰度前十的物种注释
tran = TRUE # 转化为相对丰度值
# 按照分类学门(j)合并
psdata = i %>% tax_glom(taxrank = j)# 转化丰度值
if (tran == TRUE) {psdata = psdata%>% transform_sample_counts(function(x) {x/sum(x)} )
}#--提取otu和物种注释表格
otu = otu_table(psdata)
tax = tax_table(psdata)
tax[1:3,1:7]#--按照指定的Top数量进行筛选与合并
for (i in 1:dim(tax)[1]) {if (row.names(tax)[i] %in% names(sort(rowSums(otu), decreasing = TRUE)[1:Top])) {tax[i,j] =tax[i,j]} else {tax[i,j]= "Other"}
}
tax_table(psdata)= tax##转化为表格
Taxonomies <- psdata %>% psmelt()
# head(Taxonomies)
Taxonomies$Abundance = Taxonomies$Abundance * 100
整理成facet需要的格式
这里的格式也很简单,就是需要一列“id”,这里我们将样本名修改为id,即可
# colnames(Taxonomies)[1] = "id"
Taxonomies$OTU = NULL
colnames(Taxonomies)[1] = "id"
保证颜色填充独立性
因为我们颜色填充有好几种方式,所以需要对每种颜色填充保重独立性,使用ggnewscale,这也是Y叔的包。
library(ggnewscale)
p <- p + new_scale_fill()
p
分面组合树和柱图
p3 <- facet_plot(p, panel = 'Stacked Barplot', data = Taxonomies, geom = geom_barh,mapping = aes(x = Abundance, fill = as.factor(Phylum)),color = "black",stat='identity' )
p3
修改配色
colbar <- dim(unique(dplyr::select(Taxonomies, one_of(j))))[1]
colors = colorRampPalette(c("#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD","#AD6F3B", "#673770","#D14285", "#652926", "#C84248","#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(colbar)
p3 + scale_fill_manual(values = colors)
ggtree调整布局
修改layout,设置中空等。
p = ggtree(hc,layout="fan", branch.length = "none", ladderize = FALSE) %<+% dd + geom_tippoint(size=5, shape=21, aes(fill=factor(Group), x=x)) + geom_tiplab(aes(color = Group,x=x*1.2), hjust=1)
p = p + xlim(-4,NA)
p
添加样本其他信息
如添加样品测序量柱状图、数值标签
p <- ggtree(hc) + theme_tree2()
p
head(dd)
dd$sequencenum = sample_sums(ps)
dd
data = data.frame(id = row.names(dd),sequencenum = dd$sequencenum )
head(data)
# p3
#---------添加序列
p2 <- facet_plot(p, panel = 'Number Barplot', data = dd , geom = geom_barh,mapping = aes(x = sequencenum ,fill = Group),stat='identity' )
p2
facet_plot(p2, panel='Stacked Barplot',data=dd, geom=geom_text, mapping=aes(x=sequencenum+20, label=sequencenum))
树+柱+堆叠图组合
p3 <- facet_plot(p2, panel = 'Abundance Barplot', data = Taxonomies, geom = geom_barh,mapping = aes(x = Abundance, fill = as.factor(Phylum)),color = "black",stat='identity' )
p3
撰文:五谷杂粮
责编:刘永鑫 中科院遗传发育所
猜你喜欢
- 10000+: 菌群分析
宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊 肠道指挥大脑 - 系列教程:微生物组入门 Biostar 微生物组 宏基因组
- 专业技能:生信宝典 学术图表 高分文章 不可或缺的人
- 一文读懂:宏基因组 寄生虫益处 进化树
- 必备技能:提问 搜索 Endnote
- 文献阅读 热心肠 SemanticScholar Geenmedical
- 扩增子分析:图表解读 分析流程 统计绘图
- 16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
- 在线工具:16S预测培养基 生信绘图
- 科研经验:云笔记 云协作 公众号
- 编程模板: Shell R Perl
- 生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外5000+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
学习扩增子、宏基因组科研思路和分析实战,关注“宏基因组”
点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA
ggplot2版聚类物种丰度堆叠图相关推荐
- 16S 物种丰度热图学习
### 1. 关于热图的用途(参考http://www.360doc.com/content/17/0729/17/45848444_675155815.shtml) 以RNA-seq为例,热图可以: ...
- R统计绘图-随机森林分类分析及物种丰度差异检验组合图
此文主要涉及随机森林组间变量重要性和物种丰度差异检验绘图,包含以下几部分内容: 1)随机森林分类: 2)随机森林分类变量重要性绘图: 3)物种丰度差异检验绘图 4)随机森林分类变量重要性及物种丰度差异 ...
- 物种丰度排序堆积柱形图及处理间各物种差异分析
物种丰度排序堆积柱形图及处理间各物种丰度非参数检验多组比较的R图形可视化 再美的可视化图形若缺少了统计检验就失去了灵魂而变得华而不实 测试数据及代码链接:https://pan.baidu.com/s ...
- 使用R语言获得16S物种丰度
还是获得16S物种丰度得老问题,最近在一台新机器上安装qiime1,发现有报错,对于这种停止维护的软件,也是正常现象吧,于是想别的办法解决,恰巧最近读R几本R语言的入门书,发现prop.table() ...
- R堆叠柱状图各成分连线画法:突出展示组间物种丰度变化
作者:朱微金 李陈浩 堆叠柱状图连线画法 提出问题 18年1月29日宏基因组转载了中科院生态中心邓晔组的文章<土壤细菌定量方法结合相对丰度分析揭示种群的真实变化 >.其中的图3基于堆叠柱状 ...
- 相对丰度柱状图matlab,R堆叠柱状图各成分连线画法:突出展示组间物种丰度变化...
作者:朱微金 李陈浩 堆叠柱状图连线画法提出问题 18年1月29日宏基因组转载了中科院生态中心邓晔组的文章<土壤细菌定量方法结合相对丰度分析揭示种群的真实变化 >.其中的图3基于堆叠柱状图 ...
- 按照物种丰度对OTU表格进行拆分-丰富和稀有物种识别
稀有物种 (rare taxa, RT),在所有的样本中丰度均低于0.1%: 丰富物种 (abundant taxa, AT),在所有的样本中丰度均高于1%: 中等物种 (moderate taxa, ...
- qiime2+biom+qiime1获得16S物种丰度
我们知道,不管是16S等扩增子测序,还是宏基因组,最后最重要的结果,就是物种的丰度情况了,qiime2给出的16S丰度结果是一个计数,对于许多软件来说这是可用的,那么如果我们想获得一个直接的百分比数据 ...
- R统计绘图 | 物种组成堆叠柱形图(绝对/相对丰度)
一.数据准备 数据使用的不同处理土壤样品的微生物组成数据,包含物种丰度,分类单元和样本分组数据.此数据为虚构,可用于练习,请不要作他用. # 1.1 设置工作路径 #knitr::opts_knit$ ...
最新文章
- 计算机网络学习--交换机和路由器转发数据原理
- Linux中sort、uniq、cut、wc命令详解
- Android Ion 框架 文件下载
- Vue Bootstrap OSS 实现文件上传
- 查天气43课-46课
- 数据结构(三):非线性逻辑结构-树
- pyinstaller 32位 64位的问题
- ansys16.0安装教程
- 数学建模美赛该如何准备?
- 密码库LibTomCrypt学习记录——(2.13)分组密码算法的工作模式——CCM加密认证模式
- [4G5G专题-4]:RRU 全面了解什么是4G+5G RF静态射频共享?
- 34套Java项目教程+源码包含Java swing项目 Java web项目 Java控制台项目(视频教程+源码)
- 输入某辆小轿车三次的 耗油量(升)和行驶里程(公里),计算平均油耗(升/百公里)。
- 上位机软件开发项目案例(一)_C#开发
- html实现京东网页
- 有限差分法上-椭圆系统
- 优锘科技:企业架构管理平台荣获2021年度创新产品奖
- Arduino与Proteus仿真实例-简单信号频率计数仿真
- 均线多头排列选股公式,选出均线多头刚起步的标的
- 鱼骨图和甘特图图表合集PPT模板-优页文档