論文
Reduced diversity and altered composition of the gut microbiome in individuals with myalgic encephalomyelitis/chronic fatigue syndrome
本地文件Giloteaux2016_Article_ReducedDiversityAndAlteredComp.pdf
image.png
今天的推文我們來重復(fù)一下論文中的Figure4a
數(shù)據(jù)下載鏈接 https://github.com/Nick243/Create-Giloteaux-2016-Phyloseq-Object
首先是安裝phyloseq這個(gè)包
BiocManager::install("phyloseq")
BiocManager::install("Rhdf5lib")
讀取數(shù)據(jù)
ps<-readRDS("ps_giloteaux_2016.rds")
對數(shù)據(jù)進(jìn)行預(yù)處理
這部分代碼就不介紹了厂抽，主要是為了拿到作圖數(shù)據(jù)就可以了
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這個(gè)數(shù)據(jù)集來作圖
library(ggplot2)
ggplot(data=final_df,
aes(x=variable,y=value,fill=Phylumn))+
geom_bar(stat = "identity",
position = "stack")
image.png
接下來進(jìn)行美化
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
這個(gè)對應(yīng)的是論文中對照的那個(gè)圖罕扎，這里配色不一樣择诈，因?yàn)轭伾容^多，不想在一個(gè)一個(gè)顏色單獨(dú)摘了。
最后是拼圖
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ù)據(jù)和代碼可以給推文贊賞1元獲取。贊賞了如果沒有收到回復(fù)可以加我的微信催我梨熙，我的微信是 mingyan24
歡迎大家關(guān)注我的公眾號
小明的數(shù)據(jù)分析筆記本
小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子吻贿；2串结、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)、基因組學(xué)舅列、群體遺傳學(xué)文獻(xiàn)閱讀筆記肌割；3、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記帐要！
