跟着Microbiome学作图：R语言ggplot2画堆积柱形图展示微生物门水平的相对丰度

用户7010445

发布于 2021-11-23 15:36:10

2K0

发布于 2021-11-23 15:36:10

论文

Reduced diversity and altered composition of the gut microbiome in individuals with myalgic encephalomyelitis/chronic fatigue syndrome

本地文件Giloteaux2016_Article_ReducedDiversityAndAlteredComp.pdf

今天的推文我们来重复一下论文中的 Figure4a

代码主要参考链接 https://www.nicholas-ollberding.com/post/introduction-to-the-statistical-analysis-of-microbiome-data-in-r/

数据下载链接 https://github.com/Nick243/Create-Giloteaux-2016-Phyloseq-Object

首先是安装phyloseq这个包

BiocManager::install("phyloseq")
BiocManager::install("Rhdf5lib")

读取数据

ps<-readRDS("ps_giloteaux_2016.rds")

对数据进行预处理

这部分代码就不介绍了，主要是为了拿到作图数据就可以了

ps<-readRDS("ps_giloteaux_2016.rds")
phyloseq::sample_sums(ps)
sort(phyloseq::sample_sums(ps))
(ps <- phyloseq::subset_samples(ps, phyloseq::sample_sums(ps) > 5000)) 
(ps <- phyloseq::prune_taxa(phyloseq::taxa_sums(ps) > 0, ps)) 

phyloseq::sample_data(ps)$Status <- ifelse(phyloseq::sample_data(ps)$Subject == "Patient", "Chronic Fatigue", "Control")
phyloseq::sample_data(ps)$Status <- factor(phyloseq::sample_data(ps)$Status, levels = c("Control", "Chronic Fatigue"))
ps %>% 
  sample_data %>%
  dplyr::count(Status)
table(phyloseq::tax_table(ps)[, "Phylum"])
ps_rel_abund = phyloseq::transform_sample_counts(ps, function(x){x / sum(x)})
phyloseq::otu_table(ps)[1:5, 1:5]
phyloseq::otu_table(ps_rel_abund)[1:5, 1:5]

#phyloseq::plot_bar(ps_rel_abund, fill = "Phylum")

ps_rel_abund@otu_table %>% dim()
ps_rel_abund@tax_table %>% head()
ps_rel_abund@tax_table %>% dim()
ps_rel_abund@sam_data %>% head()
ps_rel_abund@phy_tree
ps_rel_abund@refseq

ps_rel_abund@otu_table %>% class()
ps_rel_abund@otu_table %>% as.data.frame() -> df1
ps_rel_abund@tax_table %>% as.data.frame() -> df2
rownames(df2) == rownames(df1)
df1$Phylumn<-df2$Phylum
table(df1$Phylumn)

ps_rel_abund@sam_data %>% as.data.frame() -> df3
df4<-data.frame(sample_id=rownames(df3),
                sample_group=df3$Subject)
head(df4)

df1 %>% reshape2::melt(id.vars="Phylumn") %>% 
  merge(.,df4,by.x="variable",by.y="sample_id") -> final_df

接下来是用 final_df这个数据集来作图

library(ggplot2)

ggplot(data=final_df,
       aes(x=variable,y=value,fill=Phylumn))+
  geom_bar(stat = "identity",
           position = "stack")

接下来进行美化

final_df %>% 
  filter(sample_group=="Control") %>% 
  group_by(Phylumn,variable,sample_group) %>% 
  summarise(value_1=sum(value)) %>% 
  drop_na(Phylumn) -> dfa

dfa$Phylumn<-factor(dfa$Phylumn,
                    levels = names(table(dfa$Phylumn))[c(2,5,7,9,1,8,4,6,3)])  


dfa %>% 
  filter(Phylumn=="Bacteroidetes") %>% 
  arrange(value_1) -> dfa.1

dfa$variable<-factor(dfa$variable,
                       levels = rev(dfa.1$variable))

dfa %>% 
  ggplot()+
  geom_bar(aes(x=variable,y=value_1,
               fill=Phylumn),
           stat="identity",
           position = "stack")+
  scale_fill_brewer(palette = "Set1")+
  theme_minimal()+
  scale_y_continuous(expand = c(0,0))+
  theme(axis.text.x = element_blank(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line())+
  labs(x="CONTROLS",
       y="Relative Abundance (%)")

这个对应的是论文中对照的那个图，这里配色不一样，因为颜色比较多，不想在一个一个颜色单独摘了。

最后是拼图

final_df %>% 
  filter(sample_group=="Control") %>% 
  group_by(Phylumn,variable,sample_group) %>% 
  summarise(value_1=sum(value)) %>% 
  drop_na(Phylumn) -> dfa

levels<-c("Bacteroidetes","Firmicutes","Proteobacteria",
          "Verrucomicrobia",
          "Actinobacteria","Tenericutes",
          "Euryarchaeota","Fusobacteria","Cyanobacteria" )
dfa$Phylumn<-factor(dfa$Phylumn,
                    levels = levels)  



dfa %>% 
  filter(Phylumn=="Bacteroidetes") %>% 
  arrange(value_1) -> dfa.1

dfa$variable<-factor(dfa$variable,
                       levels = rev(dfa.1$variable))

dfa %>% 
  ggplot()+
  geom_bar(aes(x=variable,y=value_1,
               fill=Phylumn),
           stat="identity",
           position = "stack")+
  scale_fill_brewer(palette = "Set1")+
  theme_minimal()+
  scale_y_continuous(expand = c(0,0))+
  theme(axis.text.x = element_blank(),
        axis.line.y = element_line(),
        axis.ticks.y = element_line())+
  labs(x="CONTROLS",
       y="Relative Abundance (%)") -> pa

table(final_df$sample_group)

final_df %>% 
  filter(sample_group=="Patient") %>% 
  group_by(Phylumn,variable,sample_group) %>% 
  summarise(value_1=sum(value)) %>% 
  drop_na(Phylumn) -> dfb

dfb$Phylumn<-factor(dfb$Phylumn,
                    levels = levels)  


dfb %>% 
  filter(Phylumn=="Bacteroidetes") %>% 
  arrange(value_1) -> dfb.1


dfb$variable<-factor(dfb$variable,
                     levels = rev(dfb.1$variable))

dfb %>% 
  ggplot()+
  geom_bar(aes(x=variable,y=value_1,
               fill=Phylumn),
           stat="identity",
           position = "stack")+
  scale_fill_brewer(palette = "Set1")+
  theme_minimal()+
  scale_y_continuous(expand = c(0,0))+
  theme(axis.text = element_blank(),
        axis.line = element_blank(),
        axis.ticks = element_blank())+
  labs(x="ME/CFS",
       y=NULL) -> pb

library(patchwork)

pa+pb+plot_layout(guides = "collect")