metagenomics中，寻找区分两个类别的marker是很常见的分析。根据很多文章中的结果，摸索了一下怎么实现

数据准备

Otu为列名，添加一个分类列如status

数据

代码部分

首先我定义了一个自己常用的主题

my_theme <- function() {
    library(ggplot2)
    library(ggthemr)
    theme_set(theme_classic() +
                  theme(axis.text = element_text(size = 8, face = "bold"), #坐标轴标签大小，加粗
                        axis.title = element_text(size = 12, face = "bold"), #坐标轴标题大小，加粗
                        plot.margin = unit(rep(1, 4), "cm"), #边距，顺序：上右下左
                        panel.grid = element_blank())) #不要横向和纵向格子线
}

正儿八经开始画图

setwd("D:/bioinfo/R learn")
library(reshape2)
library(randomForest)
library(dplyr)
library(ggplot2)
library(magrittr)
library(patchwork)
my_theme()
otu <- read.table("randomforest.txt", header = T, row.names = 1, sep = "\t")
set.seed(123)
nsamp <- floor(nrow(otu))
# 随机选择70%数据作为训练集
indexes <- sample(1:nrow(otu), size = nsamp*0.7)
training <- otu[indexes,]
training$Status <- factor(training$Status)
validation <- otu[-indexes,]
validation$Status <- factor(validation$Status)
# 训练
rf_classifier <- randomForest(Status ~ ., data=training, ntree=100, importance=TRUE)
rf_classifier
varImpPlot(rf_classifier)
# 当然可以提取出来自己画图
# 首先寻找合适个数的预测变量
# 通过交叉验证来寻找
otu_train.cv <- replicate(5, rfcv(training[-ncol(training)], 
                                  training$Status, 
                                  cv.fold = 10, step = 2), simplify = FALSE)
otu_train.cv
#提取验证结果绘图
otu_train.cv.df <- data.frame(sapply(otu_train.cv, '[[', 'error.cv'))
colnames(otu_train.cv.df) <- c("err1", "err2", "err3", "err4", "err5")
otu_train.cv.df %<>% .[-nrow(.),] %>% mutate(errmean = rowMeans(.)) %>% mutate(num = as.numeric(rownames(.)))
# 选择9个预测变量
p1 <- ggplot() +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$err1), colour = 'grey', lwd = 1.5) +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$err2), colour = 'grey', lwd = 1.5) +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$err3), colour = 'grey', lwd = 1.5) +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$err4), colour = 'grey', lwd = 1.5) +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$err5), colour = 'grey', lwd = 1.5) +
    geom_line(aes(x = otu_train.cv.df$num, y = otu_train.cv.df$errmean), colour = 'black', lwd = 1.5) +
    geom_vline(xintercept = 9, colour='black', lwd=1, linetype="dashed") +
    labs(title=paste('Training set (n = ', length(indexes),')', sep = ''), 
         x='Number of classes ', y='Cross-validation error rate') + 
    coord_trans(x = "log2") +
    scale_x_continuous(breaks = c(1, 2, 5, 10, 20, 30, 50, 100, 140)) +
    annotate("text", x = 9, y = max(otu_train.cv.df$errmean), label=paste("optimal = ", 9, sep=""))
# 选择前9个important feature画图
imp <- importance(rf_classifier) %>% 
    data.frame() %>% 
    arrange(desc(MeanDecreaseAccuracy)) %>% 
    head(n = 9) %>% 
    select(3) %>%
    mutate(class = rownames(.)) %>% 
    mutate(phylum = sapply(strsplit(rownames(.), ".", fixed = TRUE), "[[",2))
imp_long <- melt(imp, id.vars = c("class", "phylum")) %>% arrange(value)
imp_long$class <- factor(imp_long$class, levels = imp_long$class)
imp_long$phylum <- factor(imp_long$phylum)

p2 <- ggplot(imp_long, aes(x= class, y = value, fill = phylum)) +
    geom_bar(stat= "identity") + 
    coord_flip() + 
    scale_fill_brewer(palette = "Paired") +
    labs(y = "Mean of decreased accuracy")
# 拼图
p1+p2+plot_annotation(title = "Fig 1", tag_levels = "A")

results

Fig 1. A 部分代码参考了如下文章，如果用到，请正确引用别人的文章（这不是我的文章）
ref: https://www.nature.com/articles/s41587-019-0104-4