跟着Nature Microbiology学作图（一）| 系统发育树绘制

一边学习，一边总结，一边分享！

引言

我们前面分享了NM的文章和代码（组学绘图代码分享，作者提供全部图的绘制代码，你还不动手收藏呢？），类似的文章和分享很多，也许你平时已经收藏了很多很多这样的教程。但是，可能你收藏后就“吃灰”了，基本没有太多时间去了解和学习。

每个人精力和时间是有限的，我们当下也许只能做好1-2件事情。因此，我们会花费很少的时间，做每期教程的分享，希望可能帮助你。

本期教程图形

获得本期教程数据和代码，在后台回复：20250207

绘图代码

加载R包

 library(ape)
library(ggnewscale)
library(ggplot2)
library(ggtree)
library(ggtreeExtra)
library(RColorBrewer)
library(tidyverse)

source("plot_functions.R")

导入数据

##'树文件
BacTree <- read.tree("pMAGs_bact_gtdtk_midroot.tree")
head(BacTree)

这是一个树形结构文件，表示系统发育树。通常，这种文件是以Newick格式或Nexus格式存储的，包含了物种或基因间的进化关系。

# 分类信息的数据表
dat <- read_tsv("pMAGS_tax.tsv")
head(dat)

这是一个包含分类信息的数据表

dataset <- read_tsv("pMAGS_Presence_Datasets.txt")
head(dataset)

covmax <- read_tsv("pMAGs_cov_sum.txt")
head(covmax)

数据处理

dat$p_c <- if_else(dat$p == "p__Proteobacteria", dat$c, dat$p)
dat$p_c <- gsub(".__", "", dat$p_c, perl = T)
dat$p_c <- gsub("_.$", "", dat$p_c, perl = T)

bacDat <- dat %>%
  filter(d == "d__Bacteria") %>%
  mutate(Abundance  = 0.2)

mostAbundantTax <- bacDat %>%
  group_by(p_c) %>%
  summarise(total   = n()) %>%
  arrange(desc(total)) %>%
  slice(1:16)

list <- mostAbundantTax$p_c
list <- c(list, "Nitrospirota")


bacDat  <- bacDat %>%
  mutate(p_c    = if_else(p_c %in% list,  p_c, "Other"))

数据整合

bacDatset <- dataset %>%
  filter(MAGs %in% bacDat$MAGs)

bacDatset$Busi <-
  gsub("Busi", "Busi et al., Nat Com, 2022", bacDatset$Busi)
bacDatset$ENSEMBLE <-
  gsub("ENSEMBLE", "Michoud et al, L&O, 2023", bacDatset$ENSEMBLE)
bacDatset$Tibet <-
  gsub("Tibet", "Tibetan Glacier Genome and Gene", bacDatset$Tibet)
bacDatset$Tara <- gsub("Tara", "Tara Oceans", bacDatset$Tara)

bacDatset <- as.data.frame(bacDatset)
rownames(bacDatset) <- bacDatset$MAGs
bacDatset$MAGs <- NULL


bactcov <- covmax %>%
  filter(MAGs %in% bacDat$MAGs) %>%
  select(MAGs, count)

bactcov <- as.data.frame(bactcov)
rownames(bactcov) <- bactcov$MAGs
bactcov$MAGs <- NULL

bactcov$count <- log10(bactcov$count)

绘图

a   <- split(bacDat$MAGs, bacDat$p_c)

tree    <- groupOTU(BacTree, a)

getPaletteBact = colorRampPalette(brewer.pal(9, "Set1"))

bactColor   <- getPaletteBact(length(unique(bacDat$p_c)) + 1)

bactColor[1]    <- "black"

p1  <-
  ggtree(tree,
         layout = 'circular',  # 表示将树形图绘制为圆形布局
         aes(color = group)) + # 表示使用group列的值来为树上的每个节点着色
  geom_tree() +
  theme_tree() +
  geom_treescale(width  = 0.1) +
  scale_color_manual(values = bactColor,
                     na.value = "transparent",
                     guide = "none") +
  #geom_text2(aes(subset=!isTip,    label=node), hjust=-.3)+
  theme(legend.position = "right")

p2 <- p1 +
  new_scale_colour() +
  new_scale_fill() +
  geom_fruit(
    data = bacDat,
    pwidth  = 0.01,
    geom = geom_bar,
    mapping = aes(y = MAGs, fill = p_c, x = 1),
    # orientation="y",
    stat = "identity",
  ) +
  scale_fill_manual(values  = bactColor[-1]) +
  labs(fill = "Taxa") +
  new_scale_colour() +
  new_scale_fill()

p3 <-
  gheatmap(
    p2,
    bacDatset,
    width = 0.2,
    offset = 0.1,
    # offset=8, width=0.6,
    colnames = FALSE,
    color = NULL
  ) +
  scale_fill_discrete(na.translate = F) +
  labs(fill = "Datasets") +
  new_scale_colour() +
  new_scale_fill()