手把手教你用OTU表绘制物种组成图

欢迎关注R语言数据分析指南

❝

本节介绍如何从最原始的OTU表开始绘制物种组成图，下面通过一个具体的例子来进行展示，各位观众老爷细细品味，希望对大家有所帮助

❞

加载R包

library(tidyverse)
library(magrittr)
library(ggh4x)
library(ggsci)

定义函数

da <- function(da) {
  (da/(da %>% rowSums())) %>% return()
}

数据清洗

❝

此处主要将绝对丰度表转化为相对丰度表，并根据注释信息将数据拆分为7类

❞

otu <- read.delim('otu.xls',row.names=1) %>%
  select_if(is.numeric) %>% t() %>% da() %>% 
  t() %>% as.data.frame() %>% rownames_to_column(var="OTU") %>% 
  left_join(.,read_tsv("otu.xls") %>% select_if(~!is.numeric(.)),
            by="OTU") %>%
  separate(taxonomy,
           into=c("domain","phylum","class","order","family","genus","species"),sep=";") %>%
mutate_at(vars(c(`domain`:`species`)),~str_split(.,"__",simplify=TRUE)[,2]) %>% 
  column_to_rownames("OTU")

数据整合

table <- otu %>% select_if(~is.numeric(.)) %>% rownames_to_column("ID")
tax <- otu %>% select_if(~!is.numeric(.)) %>% rownames_to_column("ID")

定义颜色

pal <- c("#E41A1C","#1E90FF","#FF8C00","#4DAF4A","#984EA3","#40E0D0","#FFC0CB",
         "#00BFFF","#FFDEAD","#90EE90","#EE82EE","#00FFFF","#F0A3FF", "#0075DC", 
         "#993F00","#4C005C","#2BCE48","#FFCC99","#808080","#94FFB5","#8F7C00",
         "#9DCC00","#C20088","#003380","#FFA405","#FFA8BB","#426600","#FF0010",
         "#5EF1F2","#00998F","#740AFF","#990000","#FFFF00")

定义主题

theme_niwot <- function(){
  theme_bw()+
    theme(strip.background = element_rect(fill="white",color="black"),
          panel.spacing = unit(0,"lines"),
          strip.text.x = element_text(size=8,color="black"),
          axis.text.y=element_text(size=8,color="black"),
          axis.title.y = element_text(size=10,color="black"),
          legend.key=element_blank(),
          legend.text = element_text(color="black",size=8),
          legend.spacing.x=unit(0.1,'cm'),
          legend.spacing.y=unit(0.1,'cm'),
          legend.key.width=unit(0.4,'cm'),
          legend.key.height=unit(0.4,'cm'),
          legend.background=element_blank(),
          panel.grid.major=element_blank(),panel.grid.minor=element_blank())
}

属水平丰度表清洗

❝

此处主要整合属水平的丰度表，并对其进行过滤

❞

df <- table %>% left_join(.,tax %>% select(1,genus),by="ID") %>% 
  select(-ID) %>% 
  group_by(genus) %>%
  summarise(across(where(is.numeric), ~ sum(.x, na.rm=TRUE))) %>% 
  filter(genus !="") %>% 
  column_to_rownames(var="genus") %>% 
  mutate(sum = rowSums(.), persent = sum / sum(sum) * 100) %>% 
  filter(persent > 1) %>% rownames_to_column(var="genus") %>% 
  select(-sum,-persent) %>% 
  pivot_longer(-genus) %>%
  set_colnames(c("genus","SamplesID","value")) %>% 
  right_join(.,read_tsv("group.xls",col_name=F) %>% set_colnames(c("SamplesID","group")),by="SamplesID") %>% 
  filter(value !=0)

定义因子

df$group <- factor(df$group,levels = df$group %>% as.data.frame() %>% 
                     distinct() %>% pull())

数据可视化

df %>% ggplot(.,aes(SamplesID,value,fill=genus))+
  geom_bar(stat="identity",position = "fill")+
  facet_nested(.~group,drop=T,scale="free",space="free",switch="y")+
  labs(x=NULL,y=NULL)+
  scale_y_continuous(expand=c(0,0),labels=scales::percent)+
  scale_x_discrete(expand=c(0.1,0.1)) +
  labs(fill="genus")+
  scale_fill_manual(values = pal)+
  theme_niwot()+
  theme(axis.text.x=element_text(angle = 90,color="black",vjust=0.5))

数据获取

❝

本节的内容到此结束，整篇文档从数据清洗到数据可视化可以说是细节满满，可以看到要从头绘制一张物种组成图还是很繁琐的，喜欢的小伙伴欢迎转发此文档附上一句话到朋友圈「30分钟后台截图给我」，即可获取对应的数据及代码,如未及时回复可添加我的微信

❞

欢迎大家扫描下方二位码加入「QQ交流群」,与全国各地上千位小伙伴交流

「关注下方公众号下回更新不迷路」，如需要加入微信交流群可添加小编微信，请备注单位+方向+姓名

往期推荐