其他
让maptree动起来
maptree以其优异的可视化展示方式，最近似乎博得很多人的眼球，如果我门让maptree动起来呢？
这种方式用来展示微生物群落的从属关系十分合适，下面让我们来试试：
很少的代码即可实现如此复杂的图片，比maptree简单多了。
rm(list=ls())
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"))
}
# Circlepacker package
library(circlepackeR)
# devtools::install_github("jeromefroe/circlepackeR") # If needed
# install.packages("C:/Users/liulanlan/Desktop/circlepackeR-master/", repos = NULL, type = "source")
# create a nested data frame giving the info of a nested dataset:
ps = readRDS("./ps_liu.rds")
ps1_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps1_rela
# 提取otu表格
otu_table = as.data.frame(t(vegan_otu(ps1_rela)))
otu_table$mean = rowMeans(otu_table)
otu_table$ID = row.names(otu_table)
head(otu_table)
#按照从大到小排序
otu_table<- arrange(otu_table, desc(mean))
subtab = head(otu_table,300)
dim(subtab)
# row.names(subtab) = subtab$ID
#对phyloseq取子集
ps_sub <- ps1_rela %>%
subset_taxa(
row.names(tax_table(ps1_rela ))%in%subtab$ID
)
ps_sub
#--------------------构造边和节点文件-----------------
#注意我们本次制作的是OTU水平的maptree
tax = as.data.frame(vegan_tax(ps_sub))
otu_table = as.data.frame(t(vegan_otu(ps_sub)))
head(otu_table)
#计算全部均值
otu_table$mean = rowMeans(otu_table)
#组合结果
data1 = merge(otu_table,tax,by = "row.names", all = TRUE)
head(data1)
row.names(data1) = data1$Row.names
data1$Row.names = NULL
data1$ID = row.names(data1)
data1$pathString <- paste("world", data1$Phylum, data1$Class, data1$Order, data1$Order, data1$Family, data1$Genus, data1$Species, data1$ID, sep = "/")
population1 <- as.Node(data1)
# Change the format. This use the data.tree library. This library needs a column that looks like root/group/subgroup/..., so I build it
library(data.tree)
# Make the plot
#library(RColorBrewer)#调色板调用包
#调用所有这个包中的调色板
display.brewer.all()
mi = brewer.pal(9,"Set1")
library("scales")
show_col(mi)
# You can custom the minimum and maximum value of the color range.
p <- circlepackeR(population1, size = "mean", color_min = "white", color_max = mi[1])
p