ggplot2版聚類物種丰度堆疊圖

寫在前面

隨著研究的逐漸深入，我們對繪圖的要求越來越高，各種之前使用的較少的圖形如今追求熱度和新穎程度，都開始逐漸在大文章中顯現。如下圖。

這是最近剛發表於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 & Evolution4, 122-131, doi:10.1038/s41559-019-1063-3 (2020).

這次的聚類加物種丰度展示讓我們學習一波。之前推出了用R語言的plot繪製的教程。

• R語言繪製帶聚類樹的堆疊柱形圖

但修改細節仍比較麻煩。今天更新基於ggplot2系統的教程。

加載依賴關係

這裡的ggtree需要使用19年7月以後的版本，因為這以後的版本才支持將聚類結果轉化為樹結構。

如果你的Bioconductor版本較舊，可能一直會安裝舊版ggtree。升新方法如下：

<code>## 先卸載先前的安裝控制程序
remove.packages(c(
"BiocInstaller"
, 
"BiocManager"
, 
"BiocVersion"
))
 

## 再安裝新版程序
install.packages(
"BiocManager"
)
BiocManager::install(update=
TRUE
, ask=
FALSE
)/<code>

<code>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 
return 
(as(tax,
"matrix"
))
}/<code>

導入數據

<code># 從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)))/<code>

ggtree繪製聚類樹

<code># 樣本間距離類型：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 "complete")
# 對樹分組
clus 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 =
dd[
1
:
3
,
1
:
3
] 


# ggtree繪圖 #----
p = ggtree(hc) %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/<code>

物種組成數據

<code># 指定物種組成的選項
i = ps # 指定輸入數據
j = 
"Phylum"
 # 使用門水平繪製丰度圖表
rep
 = 
6
 # 重複數量是
6
個
Top = 
10
 # 提取丰度前十的物種註釋
tran = TRUE # 轉化為相對丰度值/<code>

<code># 按照分類學門(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 % psmelt

/<code>

整理成facet需要的格式

這裡的格式也很簡單，就是需要一列“id”，這裡我們將樣本名修改為id，即可

<code>
 colnames(Taxonomies)[1] = 
"id"

Taxonomies
$OTU
 =

colnames(Taxonomies)[1] = 
"id"
/<code>

保證顏色填充獨立性

因為我們顏色填充有好幾種方式，所以需要對每種顏色填充保重獨立性，使用ggnewscale。

<code>library(ggnewscale)
p p/<code>

分面組合樹和柱圖

<code>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
/<code>

修改配色

<code>
colbar
 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)/<code>

ggtree調整佈局

修改layout，設置中空等。

<code>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
/<code>

添加樣本其他信息

如添加樣品測序量柱狀圖、數值標籤

<code>
p
 p/<code>

<code>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
/<code>

<code>facet_plot(p2, panel=
'Stacked Barplot'
,
data
=dd, geom=geom_text, mapping=aes(x=sequencenum+
20
, label=sequencenum))/<code>

樹+柱+堆疊圖組合

<code>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
/<code>

撰文：五穀雜糧

責編：劉永鑫 中科院遺傳發育所

10000+：菌群分析 寶寶與貓狗 梅毒狂想曲 提DNA發Nature Cell專刊 腸道指揮大腦

系列教程：微生物組入門 Biostar 微生物組 宏基因組

專業技能：學術圖表 高分文章 生信寶典 不可或缺的人