本節(jié)我們來繪制微生物的物種組成圖牍汹；先通過phyloseq包整合數(shù)據(jù)垛贤，之后利用MicrobiotaProcess包繪制門水平物種組成圖，最后通過ggplot2自定義繪制物種組成圖，喜歡的小伙伴可以關(guān)注我的公眾號(hào)R語言數(shù)據(jù)分析指南逸嘀，后臺(tái)回復(fù)關(guān)鍵詞物種組成圖獲取全部數(shù)據(jù)及代碼。
加載必須R包
rm(list=ls())
library(pacman)
library(magrittr)
library(reshape2)
pacman::p_load(tidyverse,phyloseq,MicrobiotaProcess,ape)
整合數(shù)據(jù)
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包進(jìn)行數(shù)據(jù)可視化
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")
導(dǎo)出門水平物種數(shù)據(jù)
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"
數(shù)據(jù)整合
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
數(shù)據(jù)可視化
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()
