一条代码完成六大微生物分类等级的和弦图绘制

写在前面

和弦图展示微生物组丰度等信息比较好看，我带大家来画一画。

实战

otupath = "./"
barpath = paste(otupath,"/circle_plot/",sep = "");otupath
dir.create(barpath)

library(phyloseq)
library(ggClusterNet)来自github我自己的R包
library(tidyverse)

for (i in 2:7) {
cir_plot(ps =ps,
Top = 10,
rank = i,
path = barpath)
}

函数

• ps phyloseq对象；

• Top：

  展示丰度最高的前几个微生物，默认十个；

• rank 按照某个分类等级进行合并微生物数据

• path 保存文件夹

cir_plot = function(ps =ps,
Top = 10,
rank = 7,
path = barpath
){
ps_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps_rela

ps_P <- ps_rela %>%
tax_glom_wt( rank = rank)
ps_P

otu_P = as.data.frame((vegan_otu(ps_P)))
head(otu_P)
tax_P = as.data.frame(vegan_tax(ps_P))

sub_design <- as.data.frame(sample_data(ps_P))
count2 = otu_P
#数据分组
iris.split <- split(count2,as.factor(sub_design$Group))
#数据分组计算平均值
iris.apply <- lapply(iris.split,function(x)colSums(x[]))
# 组合结果
iris.combine <- do.call(rbind,iris.apply)
ven2 = t(iris.combine)
Taxonomies <- ps %>%
tax_glom_wt(rank = rank) %>%
transform_sample_counts(function(x) {x/sum(x)} )%>%
psmelt() %>%
#filter(Abundance > 0.05) %>%
arrange(Phylum)
iris_groups<- group_by(Taxonomies, Phylum)
ps0_sum <- dplyr::summarise(iris_groups, mean(Abundance), sd(Abundance))
ps0_sum[is.na(ps0_sum)] <- 0
head(ps0_sum)
colnames(ps0_sum) = c("ID","mean","sd")

ps0_sum <- arrange(ps0_sum,desc(mean))
ps0_sum$mean <- ps0_sum$mean *100
ps0_sum <- as.data.frame(ps0_sum)
head(ps0_sum)
top_P = ps0_sum$ID[1:Top];top_P

### 开始进一步合并过滤
otu_P = as.data.frame(t(otu_P))
otu_tax = merge(ven2,tax_P,by = "row.names",all = F)
dim(otu_tax)
otu_tax$Phylum = as.character(otu_tax$Phylum)
otu_tax$Phylum[is.na(otu_tax$Phylum)] = "others"

i = 1
for (i in 1:length( otu_tax$Phylum)) {
if(otu_tax$Phylum [i] %in% top_P){otu_tax$Phylum [i] = otu_tax$Phylum [i]}

else if(!otu_tax$Phylum [i] %in% top_P){otu_tax$Phylum [i] = "others"}

}
otu_tax$Phylum = as.factor(otu_tax$Phylum)
head(otu_tax)

otu_mean = otu_tax[as.character(unique(sub_design$Group))]
head(otu_mean)
row.names(otu_mean) = row.names(otu_tax)
iris.split <- split(otu_mean,as.factor(otu_tax$Phylum))
#数据分组计算平均值
iris.apply <- lapply(iris.split,function(x)colSums(x[]))
# 组合结果
iris.combine <- do.call(rbind,iris.apply)
mer_otu_mean = t(iris.combine)

head(mer_otu_mean )

mer_otu_mean = t(mer_otu_mean)

library(statnet)
library(circlize)
library(RColorBrewer)#调色板调用包
mi_sam = brewer.pal(9,"Set1")
mi_tax = colorRampPalette(c( "#CBD588", "#599861", "orange","#DA5724", "#508578", "#CD9BCD",
"#AD6F3B", "#673770","#D14285", "#652926", "#C84248",
"#8569D5", "#5E738F","#D1A33D", "#8A7C64","black"))(length(row.names(mer_otu_mean)))
# library("scales")
# show_col(mi_sam )
# show_col(mi_tax)

grid.col = NULL
#这里设置样品颜色
grid.col[as.character(unique(sub_design$SampleType))] = mi_sam
#设置群落中物种水平颜色
# grid.col[colnames(mer_otu_mean)] = mi_tax
grid.col[row.names(mer_otu_mean)] = mi_tax

FileName2 <- paste(path,"/",rank.names(ps)[rank],"_cricle",".pdf", sep = "")
pdf(FileName2, width = 12, height = 8)

#gap.degree修改间隔，不同小块之间的间隔
circos.par(gap.degree = c(rep(2, nrow(mer_otu_mean)-1), 10, rep(2, ncol(mer_otu_mean)-1), 10),
start.degree = 180)
chordDiagram(mer_otu_mean,directional = F,diffHeight = 0.06,grid.col = grid.col, transparency = 0.5, annotationTrack =c("grid", "axis"),
preAllocateTracks = 2
)

circos.track(track.index = 1, panel.fun = function(x, y) {
circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index,
facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5))}, bg.border = NA) # here set bg.border to NA is important
circos.clear()
dev.off()
}