論文

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

本地文件Giloteaux2016_Article_ReducedDiversityAndAlteredComp.pdf

image.png

今天的推文我們來重復一下論文中的Figure4a

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

數(shù)據下載鏈接 https://github.com/Nick243/Create-Giloteaux-2016-Phyloseq-Object

首先是安裝phyloseq這個包

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

讀取數(shù)據

ps<-readRDS("ps_giloteaux_2016.rds")

對數(shù)據進行預處理

這部分代碼就不介紹了，主要是為了拿到作圖數(shù)據就可以了

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這個數(shù)據集來作圖

library(ggplot2)

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

image.png

接下來進行美化

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 (%)")

image.png

這個對應的是論文中對照的那個圖，這里配色不一樣，因為顏色比較多，不想在一個一個顏色單獨摘了。

最后是拼圖

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")

image.png

今天推文的示例數(shù)據和代碼可以給推文贊賞1元獲取。贊賞了如果沒有收到回復可以加我的微信催我，我的微信是 mingyan24

歡迎大家關注我的公眾號

小明的數(shù)據分析筆記本

小明的數(shù)據分析筆記本 公眾號 主要分享：1、R語言和python做數(shù)據分析和數(shù)據可視化的簡單小例子；2、園藝植物相關轉錄組學、基因組學、群體遺傳學文獻閱讀筆記；3、生物信息學入門學習資料及自己的學習筆記！