随机森林和特征可视化-(火柴图和玫瑰图)
写在前面2
今天补充几句话，这部分机器学习的完整内容也会整理到刘老师主编的、实验和分析方法的项目中，这里只是一个引子，供大家理解和学习。
写在前面1
现在机器学习，尤其是随机森林在微生物组学的应用越来越广泛，除了核心的机器学习部分外，为我们展现的往往是加工过的图片，其中最重要的图片之一就是关键变量的可视化，今天我带大家来展示一下关键变量的可视化部分。
实战
下面的ggClusterNet包内置了phyloseq对象数据，直接使用data(ps)即可导入，如果没有安装ggClusterNet，那也可以在公众号下面回复：数据，找到文件中的ps_liu.rds文件，使用readRDS来读入，即为同样效果。
R包和数据
library(phyloseq)
library(randomForest)
library(tidyverse)
# 包在github上不是很容易安装
#library(ggClusterNet)
# data(ps)
数据整理和随机僧林
这部分注意数据整理的格式，和vegan包一样的。
机器学习设置随机种子。
参数默认，调整参见帮助，或者网络上很多参数解释。最直白的评价指标：如果模型好，那就不用调整了。
map = as.data.frame(sample_data(ps))
head(map)
#-scaleing relative abundancce#----
ps_rela = transform_sample_counts(ps, function(x) x / sum(x) );ps_rela
mapping = as.data.frame(sample_data(ps_rela))
otutab = as.data.frame((vegan_otu(ps_rela)))
# Set classification info.
otutab$group = factor(mapping$Group)
# set random seed for reproducible reandomforest model
set.seed(315)
model_rf= randomForest::randomForest(group ~ ., data=otutab, importance=TRUE, proximity=TRUE)
print(model_rf)
重要变量提取和整理
这些OTU无非是提取出来重要性程度，然后加上丰度值和分类水平即可。下面的代码给大家一些参考。
#------------k可视化
###### inportant OTU picked out and plot
a=as.data.frame(round(importance(model_rf), 2))
a$id=row.names(a)
a2<- dplyr::arrange(a, desc(MeanDecreaseAccuracy))
row.names(a2)=a2$id
a3=head(a2,n=20)
taxonomy <- as.data.frame(vegan_tax(ps))
tax = taxonomy[rownames(a3),]
a3 = merge(a3,tax,by = "row.names",all = F)
row.names(a3) = a3$Row.names
a3$Row.names = NULL
OTU = vegan_otu(ps_rela)
### pice mapping
design = as.data.frame(sample_data(ps_rela))
#mean abundance by groups
iris.split <- split(as.data.frame(OTU),as.factor(design$Group))
iris.apply <- lapply(iris.split,function(x)colMeans(x,na.rm = TRUE))
norm2 <- do.call(rbind,iris.apply)%>% # combine result
t()
colnames(norm2) = paste(colnames(norm2),"mean",sep = "")
ind_fal = merge(a3,norm2,by = "row.names",all = F)
head(ind_fal)
火柴棒图形展示
p1 <- ggplot(a3, aes(x = MeanDecreaseAccuracy, y = reorder(id,MeanDecreaseAccuracy))) +
geom_point(size=6,pch=20,fill = "blue")+
geom_segment(aes(yend=id),xend=0,size=3)+
geom_label(aes(x =MeanDecreaseAccuracy*1.1, label = Phylum),size = 3)
p1
玫瑰图展示
a3<- dplyr::arrange(a3, desc(MeanDecreaseAccuracy))
a3$iid = paste(1:length(a3$id))
angle1 = 90 - 360 * ( as.numeric(a3$iid) - 0.5) /length(a3$id)
a3$id = factor(a3$id,levels = a3$id)
p2 = a3 %>%
ggplot(aes(x = factor(id), y = MeanDecreaseAccuracy ,label = Phylum)) +
geom_bar(stat = 'identity', position = 'dodge',fill = "blue") +
# scale_fill_manual(values = mi)+
geom_text(hjust = 0, angle = angle1, alpha = 1) +
coord_polar() +
ggtitle('') +
ylim(c(-1,5))+
theme_void()
p2