跟著Microbiome學(xué)作圖:R語言ggplot2畫堆積柱形圖展示微生物門水平的相對豐度

論文

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

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

數(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í)筆記帐要!

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末把敞,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子榨惠,更是在濱河造成了極大的恐慌奋早,老刑警劉巖盛霎,帶你破解...
    沈念sama閱讀 223,002評論 6 519
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異耽装,居然都是意外死亡愤炸,警方通過查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 95,357評論 3 400
  • 文/潘曉璐 我一進(jìn)店門掉奄,熙熙樓的掌柜王于貴愁眉苦臉地迎上來规个,“玉大人,你說我怎么就攤上這事姓建〉郑” “怎么了?”我有些...
    開封第一講書人閱讀 169,787評論 0 365
  • 文/不壞的土叔 我叫張陵速兔,是天一觀的道長墅拭。 經(jīng)常有香客問我,道長涣狗,這世上最難降的妖魔是什么谍婉? 我笑而不...
    開封第一講書人閱讀 60,237評論 1 300
  • 正文 為了忘掉前任,我火速辦了婚禮镀钓,結(jié)果婚禮上屡萤,老公的妹妹穿的比我還像新娘。我一直安慰自己掸宛,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 69,237評論 6 398
  • 文/花漫 我一把揭開白布招拙。 她就那樣靜靜地躺著唧瘾,像睡著了一般。 火紅的嫁衣襯著肌膚如雪别凤。 梳的紋絲不亂的頭發(fā)上饰序,一...
    開封第一講書人閱讀 52,821評論 1 314
  • 那天,我揣著相機(jī)與錄音规哪,去河邊找鬼求豫。 笑死,一個(gè)胖子當(dāng)著我的面吹牛诉稍,可吹牛的內(nèi)容都是我干的蝠嘉。 我是一名探鬼主播,決...
    沈念sama閱讀 41,236評論 3 424
  • 文/蒼蘭香墨 我猛地睜開眼杯巨,長吁一口氣:“原來是場噩夢啊……” “哼蚤告!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起服爷,我...
    開封第一講書人閱讀 40,196評論 0 277
  • 序言:老撾萬榮一對情侶失蹤杜恰,失蹤者是張志新(化名)和其女友劉穎获诈,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體心褐,經(jīng)...
    沈念sama閱讀 46,716評論 1 320
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡舔涎,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 38,794評論 3 343
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了逗爹。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片亡嫌。...
    茶點(diǎn)故事閱讀 40,928評論 1 353
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖桶至,靈堂內(nèi)的尸體忽然破棺而出昼伴,到底是詐尸還是另有隱情,我是刑警寧澤镣屹,帶...
    沈念sama閱讀 36,583評論 5 351
  • 正文 年R本政府宣布圃郊,位于F島的核電站,受9級特大地震影響女蜈,放射性物質(zhì)發(fā)生泄漏持舆。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 42,264評論 3 336
  • 文/蒙蒙 一伪窖、第九天 我趴在偏房一處隱蔽的房頂上張望逸寓。 院中可真熱鬧,春花似錦覆山、人聲如沸竹伸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,755評論 0 25
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽勋篓。三九已至,卻和暖如春魏割,著一層夾襖步出監(jiān)牢的瞬間譬嚣,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,869評論 1 274
  • 我被黑心中介騙來泰國打工钞它, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留拜银,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 49,378評論 3 379
  • 正文 我出身青樓遭垛,卻偏偏與公主長得像尼桶,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個(gè)殘疾皇子锯仪,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,937評論 2 361

推薦閱讀更多精彩內(nèi)容