写在前面堆叠柱状图做成环状，得益于Y叔的ggtreextra，因为我不仅仅要做成环状，还要添加聚类模块。所以就不太好办了。这种展示在样本很少的情况下其实不是很好，但是当样本很多的时候，尤其是上百哥，我们来发现一些规律的时候就显得非常重要了。 
library(RColorBrewer)#调色板调用包
library(ggClusterNet)
data(ps)
#---对于样本比较多的，尝试另外一种堆叠柱状图
Top = 10
dist = "bray"
cuttree = 3
hcluter_method = "complete"
colset3 <- c(brewer.pal(12,"Paired"),brewer.pal(9,"Pastel1"))# phyloseq(ps)对象标准化
ps1_rela = phyloseq::transform_sample_counts(ps, function(x) x / sum(x) )
# 导出OTU表
otu = as.data.frame(t(vegan_otu(ps1_rela)))
#计算距离矩阵
unif = phyloseq::distance(ps1_rela , method = dist)
# 聚类树，method默认为complete
hc <- stats::hclust(unif, method = hcluter_method )
# take grouping with hcluster tree
clus <- cutree(hc, cuttree )
# 提取树中分组的标签和分组编号
d = data.frame(label = names(clus),
member = factor(clus))
# eatract mapping file
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
# library(tidyverse)
# ggtree绘图 #----
p = ggtree(hc, layout='circular') %<+% dd +
geom_tippoint(size=5, shape=21, aes(fill= Group, x=x)) +
geom_tiplab(aes(color = Group,x=x * 1.2), hjust=1,offset=0.1) + xlim(-0.5,NA)
p展示丰度最高的一批微生物丰度# 按照分类学门(j)合并
psdata = tax_glom_wt(ps = ps1_rela,ranks = "Phylum")
# 转化丰度值
psdata = psdata%>% transform_sample_counts(function(x) {x/sum(x)} )
#--提取otu和物种注释表格
otu = otu_table(psdata)
tax = tax_table(psdata)
#--按照指定的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
Taxonomies$OTU = NULL
colnames(Taxonomies)[1] = "id"
head(Taxonomies)
dat2 = data.frame(id = Taxonomies$id,Abundance = Taxonomies$Abundance,Phylum = Taxonomies$Phylum)
head(dat2)
p2 <- p + new_scale_fill() + geom_fruit(
data=dat2,
geom=geom_bar,
mapping=aes(
x=Abundance,
y=id,
fill= Phylum
),
stat="identity",
width = 0.4,
orientation="y",
offset=0.7,
pwidth=3,
axis.params=list(
axis = "x",
text.angle = -45,
hjust = 0,
vjust = 0.5,
nbreak = 4
)
) +
scale_fill_manual(
name = "Number of interactions",
values=c(colset3),
guide=guide_legend(keywidth=0.5,keyheight=0.5,order=5)
) +
theme(
legend.background=element_rect(fill=NA),
legend.title=element_text(size=6.5),
legend.text=element_text(size=5),
legend.spacing.y = unit(0.02, "cm"),
legend.margin=margin(0.1, 0.9, 0.1,-0.9, unit="cm"),
legend.box.margin=margin(0.1, 0.9, 0.1, -0.9, unit="cm"),
plot.margin = unit(c(-1.2, -1.2, -1.2, 0.1),"cm")
)
p2
ggsave("./cs.pdf",p2,width = 14,height = 14)
|
