科研绘图系列：R语言微生物堆积图（stacked barplot）

专注收集和自写可发表的科研图形的数据和代码分享，该系列的数据均可从以下链接下载：

百度云盘链接: https://pan.baidu.com/s/1M4vgU1ls0tilt0oSwFbqYQ
提取码: 请关注WX公zhong号 生信学习者 后台发送 科研绘图 获取提取码

介绍

堆叠条形图是一种数据可视化图表，它通过将每个条形分割成多个部分来展示不同类别的数值。每个条形代表一个总体数据，而条形内的每个部分则代表该总体数据中不同子类别的数值。这种图表特别适合展示整体与部分的关系，以及各部分在整体中的比例。

特点：

1. 多部分条形：每个条形由多个不同颜色或图案的部分堆叠而成。

2. 直观展示比例：通过条形内各部分的长度，直观地展示各子类别在总体中的比例。

3. 易于比较：可以在同一图表中比较不同总体的子类别分布。

说明的问题：

1. 总体与部分的关系：展示一个总体数据由哪些部分组成，以及各部分的相对大小。

2. 各部分的比较：在同一图表中，可以比较不同总体中各部分的比例。

3. 数据分布：展示不同类别在总体数据中的分布情况，有助于识别主要类别和次要类别。

4. 趋势分析：通过观察不同时间点或不同条件下的堆叠条形图，可以分析数据的变化趋势。

加载R包

knitr::opts_chunk$set(warning = F, message = F)
library(tidyverse)
library(scales)
library(patchwork)

grp <- c("HTN", "Normal")
grp.col <- c("#EE2B2B", "#2D6BB4")

top.col <- c("#6C326C", "#77A2D1", "#FFD169", "#635F5F", "#D4D09A",
             "#993116", "#6798CE", "#146666", "#CE9924", "#6D659D",
             "#9F9B27", "#6D659D", "#9F9B27", "#C80b8A", "#2C3A89",
             "#C8C5C5", "#90E2BF", "#FDAB4D", "#F4F4E8", "#B054BF",
             "#FCE873", "#FFCCDB", "#AFD300", "#B089D8", "#F96E6F",
             "#AAD3ED", "#639BCE")

导入数据

数据可从以下链接下载：

• 百度云盘链接: https://pan.baidu.com/s/1M4vgU1ls0tilt0oSwFbqYQ

• 提取码: 请关注WX公zhong号 生信学习者 后台发送 科研绘图 获取提取码

phen <- read.csv("./inputdata/46-phenotype.csv")
phylum <- read.csv("./inputdata/46-phylum_count.csv") %>%
  tibble::column_to_rownames("tax")
family <- read.csv("./inputdata/46-family_count.csv") %>%
  tibble::column_to_rownames("tax")

函数

plot_StackedBar微生物单样本和分组堆积图，该函数主要执行以下步骤：

1. 输入处理：接受三个参数：x（包含样本ID和组信息的数据框），y（包含样本数据的矩阵或数据框），num（一个数字，指定显示的前N个分类学分类）。

2. 数据合并：将x和y数据合并，确保每个样本的分类学数据与相应的样本ID和组信息对应。

3. 计算均值：计算每个分类学分类的平均值，并按均值降序排列，选择前num个分类学分类。

4. 计算其他分类的均值：计算剩余分类学分类的组内和样本内平均值。

5. 生成分类学标签：将分类学名称中的点（.）替换为分号（;），并添加“Other”类别。

6. 生成组内均值数据：计算每个组内的平均值，并与前num个分类学分类关联。

7. 生成样本内均值数据：计算每个样本的平均值，并与前num个分类学分类关联。

8. 绘制组内堆叠条形图：使用ggplot2绘制组内堆叠条形图，显示每个组内不同分类学分类的相对丰度。

9. 绘制样本内堆叠条形图：使用ggplot2绘制样本内堆叠条形图，显示每个样本内不同分类学分类的相对丰度。

10. 组合图表：将组内和样本内堆叠条形图组合在一起，形成最终的可视化结果。

11. 输出结果：返回生成的堆叠条形图。

plot_StackedBar <- function(
  x,
  y,
  num = 20) {
    
  mdat <- x %>% 
    dplyr::select(SampleID, Group) %>%
    dplyr::inner_join(
      y %>% t() %>%
        data.frame() %>%
        tibble::rownames_to_column("SampleID"),
      by = "SampleID") 
  
  mdat.num.mean <-  mdat %>%
    dplyr::select(-c("SampleID", "Group")) %>%
    dplyr::summarise_each(mean) %>%
    tidyr::gather(key="tax", value="value") %>%
    dplyr::arrange(desc(value)) %>%
    dplyr::slice(c(1:num)) %>%
    dplyr::mutate(tax=as.character(tax), tax2=tax)
    
  mdat.num.mean.no <- mdat %>%
    dplyr::select(-c("SampleID", "Group")) %>%
    dplyr::summarise_each(mean) %>%
    tidyr::gather(key="tax", value="value") %>%
    dplyr::arrange(desc(value)) %>%
    dplyr::slice(-c(1:num)) %>%
    dplyr::mutate(tax=as.character(tax))
    
  if(nrow(mdat.num.mean.no) == 0){
    other.group.mean <- data.frame()
    other.individual <- data.frame()
  }else{
    other.group.mean <- mdat %>% 
      dplyr::select(-SampleID) %>%
      dplyr::group_by(Group) %>%
      dplyr::summarise_each(mean) %>%
      dplyr::ungroup() %>%
      tidyr::gather(key="tax", value="value", -Group) %>%
      dplyr::filter(tax%in%mdat.num.mean.no$tax) %>%
      dplyr::select(-tax) %>%
      dplyr::group_by(Group) %>%
      dplyr::summarise_each(sum) %>%
      dplyr::mutate(tax="Other")
    
    other.individual <- mdat %>%
      dplyr::select(-Group) %>%
      tidyr::gather(key="tax", value="value", -SampleID) %>%
      dplyr::filter(tax%in%mdat.num.mean.no$tax) %>%
      dplyr::mutate(SampleID=factor(SampleID, levels = as.character(phen$SampleID))) %>%
      dplyr::select(-tax) %>%
      dplyr::group_by(SampleID) %>%
      dplyr::summarise_each(sum) %>%
      dplyr::mutate(tax="Other")
    }
    
  taxonomy <- c(gsub("\\.", ";", mdat.num.mean$tax), "Other")
  
  mdat.group.mean <- mdat %>%
    dplyr::select(-SampleID) %>%
    dplyr::group_by(Group) %>%
    dplyr::summarise_each(mean) %>%
    tidyr::gather(key="tax", value="value", -Group) %>%
    dplyr::filter(tax%in%mdat.num.mean$tax) %>%
    dplyr::mutate(Group=factor(Group, levels = grp)) %>%
    dplyr::mutate(tax=gsub("\\.", ";", tax)) %>%
    rbind(other.group.mean) %>%
    dplyr::mutate(tax=factor(tax, levels = taxonomy))   
     
  mdat.individual <- mdat %>%
    dplyr::select(-Group) %>%
    tidyr::gather(key="tax", value="value", -SampleID) %>%
    dplyr::filter(tax%in%mdat.num.mean$tax) %>%
    dplyr::mutate(SampleID=factor(SampleID, levels = as.character(phen$SampleID))) %>%
    dplyr::mutate(tax=gsub("\\.", ";", tax)) %>%
    rbind(other.individual) %>%
    dplyr::mutate(tax=factor(tax, levels = taxonomy)) 
      
  p_group <- ggplot(mdat.group.mean, aes(x = Group, y = value, fill = tax)) + 
    geom_bar(stat= 'identity', position = 'fill',width = 0.95)+ 
    scale_fill_manual(values = top.col) +
    scale_y_continuous(labels = percent,
                       expand = c(0, 0)) + 
    labs(x = '', y = 'Relative Abundance', fill = NULL)+ 
    guides(fill = guide_legend(ncol = 1, bycol = TRUE, override.aes = list(size = 5))) +
    theme_bw()+
    theme(axis.title.y = element_text(face = 'bold',color = 'black',size = 14),
          axis.title.x = element_text(face = 'bold',color = 'black',size = 14,vjust = -1.2),
          axis.text.y = element_text(face = 'bold',color = 'black',size = 10),
          axis.text.x = element_text(face = 'bold',color = 'black',size = 12,
                                     angle = 45,vjust = 0.5), 
          panel.grid = element_blank(),
          legend.position = 'right',
          legend.key.height = unit(0.6,'cm'),
          legend.text = element_text(face = 'bold',color = 'black',size = 10)) 
  
  p_indi <- ggplot(mdat.individual, aes(x = SampleID, y = value, fill = tax)) + 
    geom_bar(stat= 'identity', position = 'fill', width = 0.95)+ 
    scale_fill_manual(values = top.col) +
    scale_y_continuous(labels = percent,
                       expand = c(0, 0)) +
    labs(x = '', y = 'Relative Abundance', fill = NULL)+ 
    guides(fill = guide_legend(ncol = 1, bycol = TRUE, override.aes = list(size = 5))) +
    theme_bw()+
    theme(axis.title.y = element_text(face = 'bold',color = 'black',size = 14),
          axis.title.x = element_text(face = 'bold',color = 'black',size = 14,
                                      vjust = -1.2),
          axis.text.y = element_text(face = 'bold',color = 'black',size = 10),
          axis.text.x = element_text(face = 'bold',color = 'black',size = 12,
                                     angle = 45,vjust = 0.5), 
          panel.grid = element_blank(),
          legend.position = 'right',
          legend.key.height = unit(0.6,'cm'),
          legend.text = element_text(face = 'bold',color = 'black',size = 10))  
  
  plot_res <- (p_indi + p_group) +  
    plot_layout(ncol = 2, widths = c(3, 1), 
    guides = "collect") & theme(legend.position='right')
  
  return(plot_res)
}

运行

• phylum

plot_StackedBar(
  x = phen,
  y = phylum,
  num = 20)

结果：左侧图展示样本内物种分布，右侧图展示组内物种分布，可以比较样本和组内的物种丰度分布情况，解析生物学上的差异。

• family

plot_StackedBar(
  x = phen,
  y = family,
  num = 20)

结果：左侧图展示样本内物种分布，右侧图展示组内物种分布，可以比较样本和组内的物种丰度分布情况，解析生物学上的差异。