一条代码得到微生物生态词云

写在前面

2020年2月19日，疫情还在继续，相信持续不了多久了。还在假期中，但还是可以持续发光发热的。今天我为大家带来词云，词云这种出图方式实际上再各个领域都有比较广发的引用，尤其是引文分析，广告等。出镜率较高。目前，鲜有再微生物生态学邻域运用到物种和功能可视化方向。直到我使用了MEGAN的可视化界面后，发现使用了这种可视化方式。java的可视化库是十分强大。甚至让我萌生了学习java的冲动。想想还是算啦，不具备这个能力。目前需要的是以最快捷的方式加快计算速度，让R语言中的许多计算过程加速一下。julia的编码风格很简单，于是就转入julia来学习。不讲废话了，代码和解释奉上。

微生物群落词云绘制

微生物有七级分类，按照每个分类等级进行词云的绘制，并将加入扩增子的主体流程，实现自动化过程。这里同样使用刘老师的示例数据，也就只需要运行一条代码，即可完成。

基础词云，用于描述不同分组就某个分类水平绘制词云

下面我从纲水平绘制微生物群落词云，这条函数会按照分组展示词云，分组列名必须为Group，这一点请大家注意，其他也没什么注意的。

# 接下来我门绘制门水品的物种词云
# j = "Phylum"
# j = "Class"
# j = "Order"
# j = "Family"
# j = "Genus"
# 微生物生态-扩增子数据词云
# 基础词云，用于描述不同分组就某个分类水平绘制词云

result = wordTax(otu = "./otutab.txt",tax = "./taxonomy.txt", map = "./metadata.tsv",j = "Class")

# 提取图片
p1 = result [[1]]
p1


高阶词云

将全部物种放到一起展示，有时候会因为物种的数据太多，词云太大，影响观看效果，所以运用ggplot分面功能，实现将词云按照分类分开展示。添加参数：facet = TRUE。

这里我添加保存图片参数，因为图形过大，所以调整limitsize = FALSE。

# facet模式，用于按照门水平绘制词云，也就是是纵坐标按照门水平分类，横坐标按照分组展示。

result = wordTax(otu = "./otutab.txt",tax = "./taxonomy.txt", map = "./metadata.tsv",j = "Family",facet = TRUE)

# 提取图片
p1 = result [[1]]

# 提取作图数据
data = result [[2]]

#
path = getwd()
FileName2 <- paste(path,"/word—facet.pdf", sep = "")

ggsave(FileName2, p1, width = 12, height =50,limitsize = FALSE)


主体函数


wordTax = function(otu = "./otutab.txt",tax = "./taxonomy.txt", map = "./metadata.tsv",j = "Phylum",facet = FALSE){

library(tidyverse)
library(ggwordcloud)
library("phyloseq")
#导入otu表格
otu = read.delim(otu,row.names = 1)
head(otu)
otu = as.matrix(otu)
str(otu)
#导入注释文件
tax = read.delim(tax,row.names = 1)
head(tax)
tax = as.matrix(tax)
# taxa_names(tax)

#导入分组文件
map = read.delim(map,row.names = 1)
head(map)

# #导入进化树
# tree = read.tree("./otus.tree")
# tree

ps <- phyloseq(otu_table(otu, taxa_are_rows=TRUE),
sample_data(map) ,
tax_table(tax)
# phy_tree(tree)

)
ps

ps_sub <- ps %>% tax_glom(taxrank = j) %>% transform_sample_counts(function(x) {x/sum(x)} )
ps_sub
otu_table = as.data.frame(t(vegan_otu(ps_sub )))

head(otu_table)

# 按照不同分组求取每组平均丰度列表
count = t(otu_table)
count2 = as.data.frame(count)
head(count)
#提取分组文件
sub_design <- as.data.frame(sample_data(ps))

#构造分组
aa = sub_design[,"Group"]
colnames(aa) = "Vengroup"
iris.split <- split(count2,as.factor(aa$Vengroup))

#数据分组计算平均值
iris.apply <- lapply(iris.split,function(x)colMeans(x[]))

# 组合结果
iris.combine <- do.call(rbind,iris.apply)
ven2 = t(iris.combine)
# sum(ven2[,1])#验证一下是否正确
ven2 = as.data.frame(ven2)
head(ven2)
ven2$ID =row.names(ven2)
# 提取物种注释文件：

vegan_tax <- function(physeq){
tax <- tax_table(physeq)

return(as(tax,"matrix"))
}
tax_table = as.data.frame(vegan_tax(ps_sub))
head(tax_table)

# match(row.names(tax_table),ven2$ID) 这是匹配的，所以直接合并

ven3 = cbind(tax_table,ven2)

#数据宽边长
library(reshape2)
plotdata<-melt(ven3,id.vars = c("ID",colnames(tax_table)),
variable.name='Group',
value.name='mean')
head(plotdata)

set.seed(42)
p = ggplot(plotdata, aes_string(label = j, size = "mean",color = "Phylum")) +
geom_text_wordcloud() +
scale_size_area(max_size = 15) +
theme_minimal() +
facet_wrap(~Group)
p

if(facet == TRUE){

p = ggplot(plotdata, aes_string(label = j, size = "mean",color = "Phylum")) +
geom_text_wordcloud() +
scale_size_area(max_size = 15) +
theme_minimal() +facet_grid(Phylum~Group)
# facet_wrap(~Group)
p
}

return(list(p,plotdata))
}