ggplot2绘制物种组成图

本节我们来绘制微生物的物种组成图；先通过phyloseq包整合数据，之后利用MicrobiotaProcess包绘制门水平物种组成图，最后通过ggplot2自定义绘制物种组成图，喜欢的小伙伴可以关注我的公众号R语言数据分析指南，后台回复关键词物种组成图获取全部数据及代码。

加载必须R包

rm(list=ls())
library(pacman)
library(magrittr)
library(reshape2)
pacman::p_load(tidyverse,phyloseq,MicrobiotaProcess,ape)

整合数据

otu_mat <- read.delim2("otu_table.tsv",header=T,
                       sep="\t",check.names = F,row.names = 1) %>%
  as.matrix()

tax_mat <- read.delim("taxa.xls",header=T,row.names = 1,
                      sep="\t",check.names = F) %>% as.matrix()
samples_df <- read.delim("group.xls",header = T,row.names = 1,
                         sep="\t",check.names = F)
tree <- read.tree("rooted_tree.tre")

OTU = otu_table(otu_mat,taxa_are_rows =T)
TAX = tax_table(tax_mat)
samples = sample_data(samples_df)

ps <- phyloseq(OTU,TAX,samples,tree)
ps

> ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 3761 taxa and 12 samples ]
sample_data() Sample Data:       [ 12 samples by 1 sample variables ]
tax_table()   Taxonomy Table:    [ 3761 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 3761 tips and 3760 internal nodes ]

MicrobiotaProcess包进行数据可视化

phytax <- get_taxadf(obj=ps, taxlevel=2)

phybar <- ggbartax(obj=phytax,facetNames="group", count=FALSE) +
  xlab(NULL) + ylab("relative abundance (%)")+
  theme(axis.text.x=element_text(face="plain",
                                 color="black",hjust=0.8,vjust=0.6,
                                 size=9, angle=90))+
  theme(strip.text.x = element_text(size=8, color="black",
                                    face="plain"))+
  theme(legend.position="right")
phybar

ggplot2自定义绘图

定义颜色

colors <-c("#E41A1C","#1E90FF","#FF8C00","#4DAF4A","#984EA3",
           "#40E0D0","#FFC0CB","#00BFFF","#FFDEAD","#90EE90",
           "#EE82EE","#00FFFF","#F0A3FF", "#0075DC", 
           "#993F00","#4C005C","#2BCE48","#FFCC99",
           "#808080","#94FFB5","#8F7C00","#9DCC00",
           "#C20088","#003380","#FFA405","#FFA8BB",
           "#426600","#FF0010","#5EF1F2","#00998F",
           "#740AFF","#990000","#FFFF00")

导出门水平物种数据

p <- phyloseq::otu_table(phytax) %>% as.data.frame()
write.table (p,file ="phylumt.xls", sep ="\t",col.names = NA)

将丰度小于1%的归类于others

computed_persent <- function(path) {
  data <- path %>%
    read.delim(check.names = FALSE, row.names = 1)
  data2 <- data %>%
    mutate(sum = rowSums(.), persent = sum / sum(sum) * 100, 
sum = NULL,) %>%
    rbind(filter(., persent < 1) %>% colSums()) %>%
    mutate(OTU_ID = c(data %>% rownames(), "others"))
  filter(data2[1:(nrow(data2) - 1),], persent > 1) %>%
    rbind(data2[nrow(data2),]) %>%
    select(ncol(.), 1:(ncol(.) - 2)) %>%
    set_rownames(seq_len(nrow(.))) %>%
    return()
}

path <- "phylumt.xls"

数据整合

a1 <- computed_persent(path) %>% melt()
a2 <- "group.xls" %>% read.delim()
a4 <- NULL

for (i in seq_len(nrow(a1))) { 
  a4[i] <- a2[which(a2[, 1] == a1[i, 2]), 2] }

a1[, 4] <- a4
a1

数据可视化

ggplot(a1,aes(variable,value,fill=OTU_ID))+
  geom_bar(stat="identity",position = "fill")+
  facet_grid(. ~ V4,scales = "free",space="free_x")+
  labs(x="",y="Proportions")+
  scale_fill_manual(values = colors)+labs(fill="")+
  theme(legend.title=element_blank())+
  scale_y_continuous(expand=c(0,0))+theme_bw()