隨著研究的逐漸深入，我們對繪圖的要求越來越高，各種之前使用的較少的圖形如今追求熱度和新穎程度，都開始逐漸在大文章中顯現。如下圖。
這是最近剛發表於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繪製的教程。
但修改細節仍比較麻煩。今天更新基於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)))
# 樣本間距離類型：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
這裡的格式也很簡單，就是需要一列「id」，這裡我們將樣本名修改為id，即可
# colnames(Taxonomies)[1] = "id"
Taxonomies$OTU = NULL
colnames(Taxonomies)[1] = "id"
因為我們顏色填充有好幾種方式，所以需要對每種顏色填充保重獨立性，使用ggnewscale。
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)
修改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+ 一線科研人員加入。參與討論，獲得專業解答，歡迎分享此文至朋友圈，並掃碼加主編好友帶你入群，務必備註「姓名-單位-研究方向-職稱/年級」。PI請明示身份，另有海內外微生物相關PI群供大佬合作交流。技術問題尋求幫助，首先閱讀《如何優雅的提問》學習解決問題思路，仍未解決群內討論，問題不私聊，幫助同行。
學習16S擴增子、宏基因組科研思路和分析實戰，關注「宏基因組」
點擊閱讀原文，跳轉最新文章目錄閱讀