本节内容属于思农数据工程产物，欢迎大家批评指正。
系统发育树实例可视化
想法来源：看到周通老师分享的关于利用itol实现发育树可视化的推文https://mp.weixin.qq.com/s/IhuiUPSw-ynp2U7xiNJxYg，加上正在学习ggtree的使用，所以想试着用ggtree来实现此图的绘制，检验一下学习成果，同时丰富一下武器库。
论文题目：Comparing with oxygen, nitrate simplifies microbial community assembly and improves function as an electron acceptor in wastewater treatment
拟复现图：Fig. 3A，如下图所示，数据链接：https://github.com/TongZhou2017/itol.toolkit/discussions/13
图片解读：由内及外进行解析。先说左图，最内圈为物种属水平系统发育树，以环形图进行展示，同时不考虑其进化距离，即各支等长；紧接着文字注释为其对应的属水平物种标签；标签文字上覆盖的色块为这些属水平物种所对应的门，具体所代表的门可参照图例示例Phylum；紧接着标签外由内至外的第一个圆环代表这些属水平物种所对应的纲，可参照图例示例Class;由内至外的第二个和第三个环实质上为柱状图，以环形进行展示，其值即柱子高低代表不同分组对应物种相对丰度乘以100后的平方根，做此处理类似于取对数等，缩小数据间的差异，但不会改变其相对大小，提升可视化的观感，蓝色代表NS组，红色代表OS组，如图例中左边的NS中的7.63即代表在左图中的所陈列的这些物种在组别NS中最高占比为7.63*7.63*%=58.22%；最外圈的环呈现的是通过相似性百分比（SIMPER）分析计算得出的对应物种在两组群落组成差异中的贡献度，值为贡献度*100后的结果，如左侧图例中的5即代表此物种对组间差异的平均贡献度为5%，贡献度差异通过颜色深浅来展示，贡献度越低颜色越浅。紧接着，右图构造与左图类似，区别就是右图无对应的物种标签同时图片整体感觉如一根根长短不一的针组合到一起，这个的主要原因感觉是因为为物种过多，同时避免图片尺寸过大，所以省略了物种标签，并且可视化时将宽度设置的很小即视觉上如一根根针一样。左图是丰富物种的可视化，右图是稀有物种的可视化。
此处采用的筛选标准为在任一样本中相对丰度超过1%的物种即为丰富物种，在任一样本均不超过1%同时在某一样本中低于0.01%即为稀有物种。关于这个定义，一直有个疑问，要不要考虑分组呢？丰富物种和稀有物种的确定是根据所有分组的所有样本，还是根据各自分组样本来确定呢？在此研究中，根据其提供附件和可视化思路，应该用的是前者，根据所有样本统一对物种进行划分类别。个人也偏向于统一定义，因为一般的物种注释都是所有样本统一进行注释的。但是看有的文献也发现有针对不同分组各自进行定义的，纠结怪就是总容易因为这些问题而纠结。回到该例图，借助ggtree，我们逐步实现左边例图的复现，同时解析下ggtree绘图步骤
#1、加载工作目录，导入所需r包，#清空环境rm(list = ls())#载入工作环境setwd("C:/系统发育树复现数据/10.1016j.envpol.2022.120243/Fig.3/Abundant-tree")#载入所需r包library(ape)library(ggtree)library(ggtreeExtra)library(ggplot2)library(extrafont)#将系统字体导入r环境中的包，可以在绘图时指定字体library(ggnewscale)#new_scale_fill函数library(dplyr)#extrafont如果第一次安装，先运行font_import()，将字体导入r，后可通过fonts()查询可使用哪些字体，#运行后如果显示已包含此字体但可视化时还提示系统无指定的字体，重启rstudio即可#2、载入数据，可视化准备，颜色设置otu.abundance <- read.table("metadata.txt",sep = "\t",header=TRUE,row.names=1)#丰度，分类等数据tree <- read.tree("abunt-tree.nwk")#系统发育树head(otu.abundance)
# 我们发现例图中标签用的是属水平注释，我们目前默认为otu名称，需要对节点标签进行替换,后面要匹配#根据发育树物种顺序，准备可视化注释文件dddd<-data.frame(id=tree$tip.label,otu.abundance[tree$tip.label,])#例图存在多个unassigned，我们加点自己的想法，对于属水平未注释出的结果我们用OTU名称替代dd$aa<-ifelse(dd$Genus=="Unassigned",dd$id,as.vector(dd$Genus))#定义可视化配色phycol <- NULLphycol[unique(dd$Phylum)[order(unique(dd$Phylum))]] <- c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99")classcol <- NULLclasscol[unique(dd$Class)[order(unique(dd$Class))]] <- c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#ff7f00","#cab2d6","#CCCCCC")#提取颜色小技巧，截图导入ps或者word，ppt等，利用取色器提取所需颜色，在设置文本颜色或填充颜色处查看其他颜色，即可查看其颜色代码NScol <- "#008AB8"OScol <- "#FE7B69"#3、准备工作进行结束，进行可视化#3.1定义基础绘图数据，%<+%,添加文件#更好的理解此函数，我们可以进行添加与不添加对比p <- ggtree(tree,layout = "fan", branch.length="none",ladderize = F,size = 0.6,open.angle =20 )head(p$data)
p<-p%<+%dd[,-which(colnames(dd)=="Class")]#避免后面添加class注释时列名冲突，不导入class列head(p$data)
#对比发现，此函数作用为，为绘图文件添加数据p#查看发育树
#3.2 添加节点标签，借助geom_tiplab()函数，设置节点标签显示内容分为我们设置的aa列p1 <- p+geom_tiplab(aes(label=aa),size=2.88,hjust=-0.01,offset=0,family="Arial",fontface=2)p1#标签过长，显示不全，我们后面还要继续加元素，所以可以暂时不管，如没有此需求，可通过xlim（）函数进行调整发育树大小，NA表示自动确定
#3.3 根据门水平注释为标签添加背景颜色，借助geom_fruit()#这里用p而不用p1是因为，绘图逻辑是图层往上添加，如果先添加标签，将会被背景颜色覆盖，因此先添加背景颜色，后添加文字p2 <- p+geom_fruit(geom = geom_tile,mapping = aes(fill = Phylum),width = 35,offset =1.5)+geom_tiplab(aes(label=aa),size=2.88,hjust=-0.01,offset=0,family="Arial",fontface=2)+scale_fill_manual(values=phycol)p2#调整背景色位置主要通过offset参数，背景大小通过width参数
#3.4 根据纲水平注释在标签外添加圈图，借助geom_cladelab()#将内部节点对纲水平进行匹配，存储于df1中tipnode <- p$data[which(p$data$isTip==TRUE),]#仅保留外部节点dd$nodeid <- tipnode$parent#提取外部节点对应的内部节点unique_tem <- dd %>% distinct(nodeid, .keep_all = TRUE)#删除重复的内部节点df1 <- data.frame(Class=unique_tem$Class,nodeid=unique_tem$nodeid)#确定借助geom_cladelab()函数添加纲水平环状图df1$label <- ""#复现图无标签，故设置显示标签为空白，可根据实际情况设定，其会显示在环图外围#这儿我们也进行一个对比，看看是否设置文本颜色，即参数textcolor所带来的区别p2+ geom_cladelab(data=df1,mapping =aes(node=nodeid,label=label,color=Class) ,extend = 0.5, barsize=3,offset = 38)+scale_color_manual(values = classcol)
p3<- p2+ geom_cladelab(data=df1,mapping =aes(node=nodeid,label=label,color=Class) ,extend = 0.5, barsize=3,offset = 38,textcolor="black")+scale_color_manual(values = classcol)
#仔细看应该可以很明显的发现右边class图例的区别，一个下方有对应的颜色字母，一个没有，这与默认显示标签有关，我们通过于aes函数外设置统一颜色解决此问题#3.5 根据NS和OS列分别添加环形丰度柱状图，借助geom_fruit()，geom_bar需设置统计方式，即stat参数p4 <- p3+geom_fruit(geom = geom_bar,mapping=aes(y=id,x=NS),pwidth=0.5,width=0.3,fill=NScol,stat="identity",offset = 1.9)+geom_fruit(geom = geom_bar,mapping=aes(y=id,x=OS),pwidth=0.5,width=0.3,fill=OScol,stat="identity",offset = 0.1)p4
#3.5 根据Dissimilarity添加最外围环形热图，借助geom_fruit()#借助new_scale_fill()函数进行新的填充设置，否则会使前后填充颜色一致，即不同区域p5 <- p4+new_scale_fill()+geom_fruit(geom = geom_tile,mapping = aes(fill = Dissimilarity),width = 4,offset =0.2)+scale_fill_gradientn(colours = c( "#E0E0E0", "#CCD5DC","#B8CAD6","#A7BECF","#92B4CD", "#7CAACC", "#649AC8","#478DC1","#1B7FBB","#0073B4","#0164AF"),breaks = seq(0,5,0.5),limits=c(0,5.2))+guides(fill = guide_colorbar(reverse = TRUE))+theme( legend.background=element_rect(fill=NA),legend.text = element_text(color = "black",size = 10),legend.title =element_text(size = 11,colour = "black",face = "bold"),legend.key.width = unit(0.65, "cm"),legend.key.height = unit(0.62, "cm"),legend.key.size = unit(6, "mm"), #调整图例刻度线大小legend.spacing = unit(0.1, "mm"), #调整不同图例之间间距)+geom_tree()+xlim(0,NA)p5
#3.6 导出pdf#由于我们前面指定了字体，所以这里导出时会发生字体消失的情况#借助showtext包中的showtext_auto()函数可解决此情况，未安装需先安装library(showtext)# 加载并使用 showtextshowtext_auto()="plot.pdf",p5,height=25/2.54,width=28/2.54)
到这里，我们就成功完成了fig3a左图的复现，虽然图形看上去很复杂酷炫，但是利用ggtree和ggtreeExtra进行复现的过程中，发现运行代码简单直接。这时候，不得不讲一句，Y叔yyds。在整个做的过程中发现，需要花一定的时间上手，理清逻辑后，使用起来还是非常便捷的，因此，私认为ggtree为发育树个性可视化的不二之选。图3b的复现我们就留到下期的推文吧，实际实现过程跟此类似，大家可以自己先运行此代码进行左图的学习复现，后趁热打铁实现右图的复现，不要忘了在成功做出来之后给小编点个赞哦。本文基本实现了95%的复现，剩下的在ai里添加优化一下，基本就跟原图一模一样啦，好看美观又不是很难，还不赶快跟着小编一起学起来。
作者：思农生信团队-汪志祥
