image
写在前面
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
image
高阶词云
将全部物种放到一起展示，有时候会因为物种的数据太多，词云太大，影响观看效果，所以运用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)
image
主体函数
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))}
以上就全部内容，如果运行出现问题可以在后台回复：词云，按照指示得到原始数据和测试代码。
image
网友评论