论文
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这个包
代码语言：javascript复制BiocManager::install("phyloseq")
BiocManager::install("Rhdf5lib")
读取数据
代码语言：javascript复制ps<-readRDS("ps_giloteaux_2016.rds")
对数据进行预处理
这部分代码就不介绍了，主要是为了拿到作图数据就可以了
代码语言：javascript复制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这个数据集来作图
代码语言：javascript复制library(ggplot2)
ggplot(data=final_df,
aes(x=variable,y=value,fill=Phylumn))
geom_bar(stat = "identity",
position = "stack")
接下来进行美化
代码语言：javascript复制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 (%)")
这个对应的是论文中对照的那个图，这里配色不一样，因为颜色比较多，不想在一个一个颜色单独摘了。
最后是拼图
代码语言：javascript复制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")
