stackbar.jpg

今天要學(xué)習(xí)的圖來(lái)自2021年10月29號(hào)發(fā)表在的Nature Communication上的一篇文章，題目是【新冠肺炎患者呼吸道菌群組成及其與宿主相互作用的臨床研究】。今天先來(lái)復(fù)現(xiàn)其中的一幅物種豐度堆積柱狀圖。

Snipaste_2021-11-11_00-04-49

DOI:10.1038/s41467-021-26500-8

[TOC]

22

讀圖

Snipaste_2021-11-15_21-09-21

該圖為分組的物種屬水平的相對(duì)豐度堆積柱狀圖。其中只顯示了前15豐度的屬，剩下的屬都?xì)w于others里

示例數(shù)據(jù)及作圖前準(zhǔn)備

由于作者給開(kāi)源的數(shù)據(jù)設(shè)置了權(quán)限，就用手上的Phylum水平的絕對(duì)豐度矩陣來(lái)作為示例。

導(dǎo)入數(shù)據(jù)

# 導(dǎo)入數(shù)據(jù)并查看數(shù)據(jù)集格式
rm(list = ls())
setwd("F:\\~\\mzbj\\mzbj_note\\NC\\3.stackbar")
phylum <- read.table('count_2Phylum.txt', sep="\t", header=T, row.names=1)
head(phylum)

> head(phylum)
                             KO1   KO2   KO3   KO4  KO5  KO6  OE1   OE2   OE3  OE4
Acidobacteria                 52    69    39    45   50   59  196    90   140   55
Actinobacteria              7771 13581 10242 10116 6305 8058 8130 10355 10470 8324
Armatimonadetes                2     1     1     2    0    4    2     7     2    5
Bacteroidetes                845   821  2766  1059 1889  782  975  2408  2783 2250
Candidatus_Saccharibacteria    1     1     6     0   33    4    0    10     1    0
Chlamydiae                     4     2     0     0    2    9   10     9    11   16

數(shù)據(jù)處理

# 將絕對(duì)豐度轉(zhuǎn)化為百分比形式的相對(duì)豐度
phylum_per <- as.data.frame(lapply(phylum, function(x) x / sum(x)))
row.names(phylum_per) <- row.names(phylum) #加一下行名
# 計(jì)算每個(gè)門(mén)水平的平均豐度 便于后續(xù)篩選                                 
phylum.ave <- apply(phylum_per, 1, FUN=mean)
phylum.2 <- cbind(phylum_per, phylum.ave)[order(-phylum.ave),] #排個(gè)序

# 選擇豐度最高的10個(gè)門(mén) 剩下的放入others里
phylum.2 <- subset(phylum.2[1:10,], select=-phylum.ave)
# 統(tǒng)計(jì)others豐度
phylum.2 <- rbind(phylum.2, others=apply(phylum.2, 2, function(x){1-sum(x)}))
# 加一列行名 便于后續(xù)的長(zhǎng)寬轉(zhuǎn)換
phylum.2 <- cbind(PhylumID=row.names(phylum.2), phylum.2)

# 長(zhǎng)寬轉(zhuǎn)換
library(reshape2)
# 因子排個(gè)序
phylum.2$PhylumID <- factor(phylum.2$PhylumID, levels = rev(phylum.2$PhylumID))
phylum.gg <- melt(phylum.2, id.vars="PhylumID", variable.name="SampleID", value.name="Abundance")
head(phylum.gg)

> head(phylum.gg)
        PhylumID SampleID  Abundance
1 Proteobacteria      KO1 0.66311258
2 Actinobacteria      KO1 0.25731788
3  Bacteroidetes      KO1 0.02798013
4     Firmicutes      KO1 0.01394040
5    Chloroflexi      KO1 0.01480132
6     Unassigned      KO1 0.01864238

# 添加分組信息
phylum.gg$group <- c(rep('KO', 66), rep('OE', 66), rep('WT', 66)) # 根據(jù)樣本情況設(shè)置

繪制

# 為了復(fù)現(xiàn)文章中的圖需要的顏色包
library(wesanderson)
library(colortools)
library(ggpubr)

ggbarplot(phylum.gg, x = "SampleID", y="Abundance", color="black", fill="PhylumID",
          legend="right", 
          legend.title="Top Phylum", main="Relative counts per Phylum",
          font.main = c(14,"bold", "black"), font.x = c(12, "bold"), 
          font.y=c(12,"bold")) + 
  theme_bw() +
  rotate_x_text() + 
  scale_fill_manual(values=c("gray",rev(wheel("skyblue3")[1:10]))) + # 顏色設(shè)置
  facet_grid(~ group, scales = "free_x", space='free') + 
  labs(x = "Sample", y = "Relative proportion") + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title = element_text(face = "bold"), 
        plot.title = element_text(face = "bold"), 
        legend.title = element_text(face = "bold")) 
ggsave(filename = "relative_counts.pdf", device="pdf", width=8, height=4)

結(jié)果展示

Snipaste_2021-11-15_22-52-59

數(shù)據(jù)及代碼

見(jiàn)原文：https://mp.weixin.qq.com/s/5mxpOCOnmO2prRkD6Y0Ygw

后記

關(guān)于更<u>詳細(xì)的代碼講解、作者的原代碼的一些細(xì)節(jié)以及我修改的地方</u>會(huì)在之后的視頻教程中詳細(xì)講到，有興趣的可以關(guān)注我的B站【木舟筆記】。

往期內(nèi)容

1. 跟著Nature學(xué)作圖 | 配對(duì)啞鈴圖+分組擬合曲線+分類(lèi)變量熱圖
2. (免費(fèi)教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
3. 跟著Nat Commun學(xué)作圖 | 1.批量箱線圖+散點(diǎn)+差異分析
4. 跟著Nat Commun學(xué)作圖 | 2.時(shí)間線圖