GraPhlAn：最美进化树或层级分类树学习笔记

Graphlan可视化进化树

撰文：文涛 南京农业的大学

责编：刘永鑫 中科院遗传发育所

GraPhlAn简介

总体来讲，Graphlan是一个可视化进化树和基于分类等级绘制层级分类树的工具。可以制作分类树是它不同于R包ggtree和iTOL的地方。然就进化树而言，iTOL功能最为全面，ggtree最容易上手，从功能上来讲graphlan很多地方不如iTOL，比如，iTOL无限制添加的数据集，外环可以制作各种图形包括箱线图等，可以使用多种符号填充外环，但是Graphlan就不能这么随便了，只能使用两个符号填充外环。再多的环属性也就是设置环数量和颜色，透明度了。ggtree最容易上手，但是就一张圈图来说，它不能添加除了热图以外的其他图形，但是在非圈图的模式下，可以对多种数据进行合并，方法将更为简单，操作也容易一些。

素颜vs美化

下面通过我们的实战，将这个原始的图形通过graphlan大法，改变到下面的样子

完全认不出来了吧，这就是美颜的重要性。

认识软件参数

参数虽多，你也不需要记住，只要认识即可，而且都是字面意思，懂英文看一遍全懂了，需要时可以来查。大多时间我们是在别人基础上修改，能看懂代码即可。

熟悉grahlan注释参数-全局参数，

一般全局参数使用方式为：参数名 + 参数值

分支控制，一般不做控制 branch_bracket_depth 0##控制一个分支块的聚集程度，值越大，树枝越集中

branch_bracket_width 0 0###控制一个分支点延伸出来的分支线的长度，设定大了会反向延长超多分支点。

annotation_background_alpha ：只能设置全部的阴影的透明度 组合命令控制注释阴影的半径

一下两个指标设置为负值，且两者均设置为-0.7，会大大减小注阴影半径

annotation_background_separation -0.40

annotation_background_offset -0.34 一般全局参数 title：设置图片题目

title_font_size：60 题目字体大小

total_plotted_degrees 350：这个设置整个圆圈的展看角度，不会设置超过360度，一般设置一个比360度小一点的值是为了添加标签使用（例如：对环添加标签ring_label）。

start_rotation 270：对圆圈图进行旋转的起始位置；注意图形下方为270°，左方为0°。

class_legend_font_size 图例大小设置

## 全局注释

annotation_background_alpha 0.02：全部注释背景透明度；

annotation_background_width 0.03 注释背景宽度

annotation_font_size：设置注释字体大小

annotation_legend_font_size 设置图上容不下时，代表字母的对应注释的字体大小

branch_thickness：全局树枝的粗细设置

branch_color ：设置全局树枝颜色

ignore_branch_len 10：控制全局树伸展和大小

clade_marker_size：全局的点的尺寸

clade_marker_edge_width：全局的点的边缘宽度

clade_separation 0.5 进化树分枝簇的紧凑程度

ring_internal_separator_thickness 0.01 设置环与环之间的厚度

ring_separator_color ：设置环间隔的颜色，结合厚度设置选项增多

ring_width：设置在这个环上每个点所占有的宽度；

ring_height：设置环与环之间的距离，默认环与环之间的距离为1；

ring_shape：对分类树外面每层的空间设置形状填充，有不多的中形状可以选择。v：正v字型；^：v字型；

ring_alpha 0.01：这是对每个叶节点占据的位置不透明度，最大值为1。值越高颜色越深，当然需要结合设置颜色ring_color

ring_color：每条环的颜色

设置全局的环参数，必须建立在叶节点对应的shape上，不可只设置全局参数，尤其是全局ring color，仅仅设置该参数是无法进行全环上色的。默认节点shape是R，也就是整个叶节点对应的方块。

节点基本参数

节点参数使用方式：节点名称 + 参数名称 + 参数值

#节点基本参数

clade_marker_size 节点的大小，从1到200都可以设置

clade_marker_color：颜色：默认无色，r：红色；b：蓝色（基本为颜色首字母）

clade_marker_shape：形状：默认圆形，*：代表五角星，h：代表六边形

clade_marker_label：设置点内部写一个标签，后面跟上标签内容

clade_marker_font_color：设置点内标签的字体颜色

clade_marker_edge_color #555555 控制节点的颜色

clade_marker_edge_width 1.2#节点的边的宽度 节点对应的形状

这些标签有我们熟悉的点，圈，各个方向的三角，五边形，五角星，钻石等。具体符号和对应名称如下：

'.' : point marker

',' : pixel marker

'o' : circle marker

v' : triangle_down marker

'^' : triangle_up marker

'<' : triangle_left marker

'>' : triangle_right marker

'1' : tri_down marker

'2' : tri_up marker

'3' : tri_left marker

'4' : tri_right marker

's' : square marker

'p' : pentagon marker

'*' : star marker

'h' : hexagon1 marker

'H' : hexagon2 marker

'+' : plus marker

'x' : x marker

'D' : diamond marker

'd' : thin_diamond marker

'|' : vline marker

'_' : hline marker 节点注释 #节点注释

annotation_background_color

annotation_background_alpha 0.02：注释背景透明度；

annotation_background_width 0.03 注释背景宽度

annotation_font_size：设置注释字体大小

annotation_legend_font_size设置图上容不下，设置成代表字母的对应注释的注释字体大小 枝条参数 branch_thickness 枝条的厚度

branch_color 枝条颜色

branch_color_from_ancestor：使用上级分支的颜色 环参数 ring_internal_separator_thickness 0.01 设置环与环之间的厚度

ring_separator_color ：设置环间隔的颜色，结合厚度设置选项增多

ring_width：设置在这个环上每个点所占有的宽度；

ring_heigh：设置环与环之间的距离，默认环与环之间的距离为1；

ring_shape：对分类树外面每层的空间设置形状填充，有不多的形状可以选择。v：正v字型；^：^字型；

ring_alpha 0.01：这是对每个叶节点占据的位置不透明度，最大值为1。值越高颜色越深，当然需要结合设置颜色ring_color

ring_color：每条环的颜色 进化树实战

使用材料

1. 我们使用已经做好的16s序列构建的进化树（由Fastree生成，一般为rep_set.tree）
2. 对应OTU的丰度表格和注释文件，分别为otu_table.txt和rep_seqs_tax.txt，有时可能两者存于同一文件内；
3. 实验设计，即样本对应属性，如分组信息(可选)。

分析要求

1. 做一张进化树，并将不同门的OTU节点标注不同的颜色和形状
2. 对不同的门的OTU添加不同颜色的阴影(低透明度颜值高)。
3. 根据每个样品的OTU丰度使用不同颜色，添加外环。
4. 全部样品分为两个组，根据OTU在两个组中的丰度信息添加不同的标签环
5. 进行差异分析，添加外环对差异OTU根据其分组进行上色并标注。
6. 将OTU总体丰度通过最外环高度展示出来。

上面我们见过软件的参数，每个节点都要注释，如果手动搞，可能真要3天至半个月画个美颜树。只要是批量操作且有规律，我们就可以用编程来实现。

使用R语言构造注释文件，这里我们明确，基于GraPhlAn的注释文件按照注释占有列数量我分为一下三个部分：

1. 占两列的注释：题目，全局参数，例如，全局点和枝参数
2. 占三列注释：叶，枝参数，全局环参数。图例参数
3. 占四列参数：环参数

由于注释文件编写耗费很大时间，Excel适合编写的也仅仅只有占两列的数据，和总环的参数。这部分编写完成之后每次绘图基本只根据不同的实验，会小范围修改。参考模板：annot1.txt

首先根据之前要求，准备要处理文件

并确定工作目录包括所需的OTU表、物种注释、实验设计等文件，文件名与代码中不对应请自行修改。

测试数据和代码下载链接，后台回复“graphlan”获取

准备需要展示的数据

OTU表和实验设计的读取、交叉筛选和标准化

# 读取OTU表

otu = read.table("otu_table.txt", sep="t",row.names= 1,header=T, check.names=F)

head(otu)

# 读入实验设计

design = read.table("metadata.txt", header=T, row.names= 1, sep="t")

head(design, n=12L)

# 交叉筛选，保持实验设计和OTU中样本对应

idx = rownames(design) %in% colnames(otu)

sub_design = design[idx,]

# 过滤OTU表

OTU = otu[,rownames(design)]

head(OTU)

OTU = as.matrix(OTU)

# 原始reads count表标准化

norm = t(t(OTU)/colSums(OTU,na=T)) #* 100 可 normalization to total 100

# 转置并转换为数据框

library(dplyr)

norm1 = norm %>% t %>% as.data.frame

norm1[1:3,1:3]

计算分组均值

# 数据分组计算平均值

otutab.split = split(norm1,as.factor(sub_design$SampleType_RS))

otutab.apply = lapply(otutab.split,function(x)colSums(x))

norm2 = do.call(rbind,otutab.apply) %>% t

# OTU表后面添加分组均值

OTU_all= cbind(as.data.frame(norm),norm2)

head(OTU_all)

追加物种注释

# 读取OTU对应的物种注释

# 9列，分别为ID，界门纲目科属种和置信度

tax = read.delim("rep_seqs_tax.txt", sep="t",row.names= 1,header=F, check.names=F)

# 只保留物种7级分类

tax = tax[1:7]

head(tax)

# 添加分类级列名：界门纲目科属种

colnames(tax) =c("kingdom","phylum","class","order","family","genus","species")

# OTU表再添加物种注释

index = merge(OTU_all,tax, by="row.names",all=F)

row.names(index) = index$Row.names

index$Row.names = NULL

head(index)

上面，我们制作OTU对应丰度、平均丰度和物种注释的数据

映射丰度至颜色

# 设定颜色梯度，橙-白-绿三色过滤

colorg = colorRampPalette(c( "#D95F02", "white","#1B9E77"))(12)

# 通用绘图缩放方法 Generic plot scaling methods

library("scales")

# 显示颜色和对应16进制RGB代码

show_col(colorg)

# 建立12种颜色的数据框

c = data.frame(id = c(1:12),col = colorg)

# 提取前面OTU中12个样品的丰度部分

a = index[1:12]

for(i in 1:nrow(a)){

a[i,1:12] = match(colnames(a[i,1:12]), colnames(sort(a[i,1:12])))

}

# 填充颜色:替换丰度顺序为对应的颜色

for(i in 1:nrow(a)){

aa = a[i,]

aa = as.data.frame(aa)

colnames(aa) = aa[1,]

ccc =t(c[colnames(aa),])

a[i,] = ccc[2,]

}

head(a)

out = cbind(index,a)

head(out)

映射分组为开关

##组比对映射形状

HH = c(rep("A",nrow(out)))

LL = c(rep("A",nrow(out)))

# 两组比较，大的标为v，小的为空

for(i in 1:nrow(out)){

if(out[i,13] > out[i,14]){

HH[i] ="v"

LL[i] =" "

}else if(out[i,13] == out[i,14]){

HH[i] =" "

LL[i] =" "

}else if(out[i,13] < out[i,14]){

HH[i] =" "

LL[i] ="^"

}

}

# 继续追加至表末

out2 = cbind(out,HH,LL)

head(out2)

差异分析

###下一步添加差异分析结果

# 安装bioconductor包安装工具

# if (!requireNamespace("BiocManager", quietly = TRUE))

# install.packages("BiocManager")

# 安装差异比较包DESeq2

# BiocManager::install("DESeq2", version = "3.8")

# BiocManager::install("GenomeInfoDb", version = "3.8")

# BiocManager::install("S4Vctors", version = "3.8")

library(DESeq2)

library(limma)

# BiocManager::install("pasilla", version = "3.8")

library(pasilla)

# library(DESeq)

head(OTU )

count = OTU

count=as.matrix(count)

#######首先比较发病和健康中的差异#############

# 读入矩阵和实验设计，指定分类

dds = DESeqDataSetFromMatrix(countData = count,

colData = sub_design,

design = ~ SampleType_RS)

dds2 = DESeq(dds)

##第二步,标准化

resultsNames(dds2)

# 将结果用results函数来获取，赋值给res变量

res = results(dds2, contrast=c("SampleType_RS","BL", "BH"),alpha=0.05)

# summary一下，看一下结果的概要信息

summary(res)

# 按校正P值和差异倍数，添加颜色

res$level = as.factor(ifelse(res$padj < 0.05 & res$log2FoldChange > 1, "#D95F02",ifelse(res$padj < 0.05 & res$log2FoldChange < -1, "#1B9E77"," ")))

head(res)

# 提取筛选的显著差异，添加至otu表

res1 = as.data.frame(res[ncol(res)])

OTU_TAX_all = merge(out2,res1, by="row.names",all=F)

head(OTU_TAX_all)

dim(OTU_TAX_all)

row.names(OTU_TAX_all) = OTU_TAX_all$Row.names

OTU_TAX_all$Row.names = NULL 构建三列注释文件，主要是节点注释 ##构建进化树节点属性注释文件

# head(OTU_TAX_all)

pcol = c(rep("A",nrow(OTU_TAX_all)))

psha = c(rep("A",nrow(OTU_TAX_all)))

# 按门水平手动设置颜色和形状

for(i in 1:nrow(OTU_TAX_all)){

if(OTU_TAX_all[i,16] == "p__Acidobacteria"){

pcol[i] ="#B0171F"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Proteobacteria"){

pcol[i] ="#D15FEE"

psha[i] ="o"

}else if(OTU_TAX_all[i,16] == "p__Acidobacteria"){

pcol[i] ="#B0171F"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Verrucomicrobia"){

pcol[i] ="#00CD00"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Cyanobacteria"){

pcol[i] ="#87CEFA"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Actinobacteria"){

pcol[i] ="#FEC80A"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Bacteroidetes"){

pcol[i] ="#EE6A50"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Chloroflexi"){

pcol[i] ="#7B68EE"

psha[i] ="D"

}else if(OTU_TAX_all[i,16] == "p__Gemmatimonadetes"){

pcol[i] ="#9ACD32"

psha[i] ="*"

}else if(OTU_TAX_all[i,16] == "p__Firmicutes"){

pcol[i] ="#8DEEEE"

psha[i] ="*"

}else {

pcol[i] ="#006400"

psha[i] ="*"

}

}

# OTUID对应物种颜色形状数据框

ano4 = cbind(row.names(OTU_TAX_all),pcol,psha)

# head(ano4)

ano4 = as.data.frame(ano4)

生成三列的graphlan要求结点注释文件

# 节点颜色映射

ste1 = data.frame(a = ano4$V1,b = c(rep("clade_marker_color",nrow(OTU_TAX_all))),c = ano4$pcol)

head(ste1)

# 节点形状映射

ste2 = data.frame(a = ano4$V1,b = c(rep("clade_marker_shape",nrow(OTU_TAX_all))),c = ano4$psha)

head(ste2)

# 节点阴影映射颜色，同节点颜色

ste3 = data.frame(a = ano4$V1,b = c(rep("annotation_background_color",nrow(OTU_TAX_all))),c = ano4$pcol)

head(ste3)

# 注释透明度，淡色，如0.1透明

ste4 = data.frame(a = ano4$V1,b = c(rep("annotation_background_alpha",nrow(OTU_TAX_all))),c = c(rep(0.1,nrow(OTU_TAX_all))))

head(ste4)

ste4$c = as.factor(ste4$c)

# 节点大小，这里统一写10，LEfSe中按丰度大小变化

ste5 = data.frame(a = ano4$V1,b = c(rep("clade_marker_size",nrow(OTU_TAX_all))),c = c(rep(10,nrow(OTU_TAX_all))))

head(ste5)

ste5$c = as.factor(ste5$c)

# 保存文件

ste_all = rbind(ste1,ste2,ste3,ste4,ste5)

head(ste_all)

# 保存统计结果，有waring正常

write.table(ste_all, "annot2.txt", append = F, quote = F, sep="t", eol = "n", na = "NA", dec = ".", row.names = F, col.names = F) 构建四列注释文件，主要是环注释 # 构建单样品丰度映射文件，共12列

# 提取之前按丰度构建的颜色矩阵

anno1 = OTU_TAX_all[22:33]

head(anno1)

# 设置样本对应环编号

colnames(anno1 ) = c(1:3,7:12,4:6)

anno1$id = row.names(anno1)

library("reshape2")

# 构建3列文件，ID，环编号，对应颜色

anno11 = melt(anno1,id.vars = "id",variable.name = "ring",value.name = "color")

# 添加第4列，类型

anno11$ringcolor = c(rep("ring_color",nrow(anno11)))

# 调整列的序列，前两列为ID和类型，其它在后面

anno11 = select(anno11,id,ringcolor,everything)

head(anno11) ##构建平均丰度标志映射文件

anno2 = data.frame(OTU_TAX_all$LL,OTU_TAX_all$HH,row.names = row.names(OTU_TAX_all))

head(anno2)

# 位置为13/14圈

colnames(anno2 ) = c(13,14)

anno2$id = row.names(anno2)

anno22 = melt(anno2,id.vars = "id",variable.name = "ring",value.name = "shape")

anno22$ringshape = c(rep("ring_shape",nrow(anno22)))

anno22 = select(anno22,id,ringshape,everything)

# 只保留有差异，即有形状标记的行

anno222 = anno22[anno22$shape %in% c("^","v"),]

head(anno222) ##构建差异分析颜色标记文件

anno3 = data.frame(id = row.names(OTU_TAX_all),

ringcolor = c(rep("ring_color",nrow(OTU_TAX_all))),

ring =c(rep(15,nrow(OTU_TAX_all))),color = OTU_TAX_all$level )

anno3$ring = as.factor(anno3$ring)

head(anno3)

anno33 = anno3[anno3$color %in% c("#D95F02","#1B9E77"),] # 合并1，3批的颜色注释

head(anno11);dim(anno11)

ann_all = rbind(anno11,anno33)

# 第二批为形状，需改列名再合并

colnames(anno222) = colnames(ann_all)

ann_all2 = rbind(ann_all,anno222)

head(ann_all2) # 添加物种丰度总体特征

# 计算所有样品的均值

norm$mean = apply(norm,1,mean)

head(norm)

# 计算平方根，标准化数据更均匀

norm$mean = sqrt(norm$mean)

norm$mean = sqrt(norm$mean)

# 第16列行高度展示丰度的4次方根

ste6 = data.frame(a = ano4$V1,b = c(rep("ring_height",nrow(OTU_TAX_all))),c = c(rep(16,nrow(OTU_TAX_all))),d = norm$mean)

head(ste6)

ste6$c = as.factor(ste6$c) # 统一列名，并合并

colnames(ste6) = colnames(ann_all2)

ann_all3 = rbind(ann_all2,ste6)

head(ann_all3)

# 保存文件

write.table(ann_all3, "annot3.txt", append = F, quote = F, sep="t", eol = "n", na = "NA", dec = ".", row.names = F, col.names = F) 总体环修饰 指标可以保存基本小修即可 title Metagenomic

title_font_size 12

total_plotted_degrees 350

start_rotation 270

branch_bracket_depth 0

branch_bracket_width 0

branch_color g

branch_thickness 0.5

ignore_branch_len 0

clade_marker_size 0.6

clade_marker_edge_width 0.35

class_legend_font_size 9

annotation_background_separation -0.4

annotation_background_offset -0.33

annotation_background_width 0.03

annotation_background_alpha 0.01

annotation_font_size 8

annotation_legend_font_size 9

ring_internal_separator_thickness 1 0.3

ring_internal_separator_thickness 2 0.3

ring_internal_separator_thickness 3 0.1

ring_internal_separator_thickness 4 0.1

ring_internal_separator_thickness 5 0.1

ring_internal_separator_thickness 6 0.1

ring_internal_separator_thickness 7 0.1

ring_internal_separator_thickness 8 0.1

ring_internal_separator_thickness 9 0.1

ring_internal_separator_thickness 10 0.1

ring_internal_separator_thickness 11 0.1

ring_internal_separator_thickness 12 0.1

ring_internal_separator_thickness 13 0.1

ring_internal_separator_thickness 14 0.1

ring_internal_separator_thickness 15 0.1

ring_internal_separator_thickness 16 0.1

ring_separator_color 1 y

ring_separator_color 2 y

ring_separator_color 3 b

ring_separator_color 4 b

ring_separator_color 5 b

ring_separator_color 6 b

ring_separator_color 7 b

ring_separator_color 8 b

ring_separator_color 9 b

ring_separator_color 10 b

ring_separator_color 11 b

ring_separator_color 12 b

ring_separator_color 13 r

ring_separator_color 14 b

ring_separator_color 15 r

ring_separator_color 16 r

ring_width 1 0.3

ring_height 1 0.5

ring_height 2 0.1

ring_height 3 0.1

ring_height 4 0.1

ring_height 5 0.1

ring_height 6 0.1

ring_height 7 0.1

ring_height 8 0.1

ring_height 9 0.1

ring_height 10 0.1

ring_height 11 0.1

ring_height 12 0.1

ring_height 13 0.1

ring_height 14 0.1

ring_height 15 0.1

ring_height 16 0.1 graphlan出图，三条命令三张图

注：如果没有graphlan，需要在linux服务器上，按官网或conda安装

推荐conda一条命令安装conda install graphlan

基于结点颜色、形状、背景色+透明度(annon2.txt)美化进化树（rep_set.tree）

# 进化树结合注释

graphlan_annotate.py --annot annot2.txt rep_set.tree 1.xml

# 绘制150dpi的PNG

graphlan.py --dpi 150 1.xml 1.png

# 绘制矢量图PDF

graphlan.py 1.xml 1.pdf

在上图的基础上，外圈添加丰度(anno3.txt)，包括12个样品，各组中高/低，极是否显著

graphlan_annotate.py --annot annot3.txt 1.xml 2.xml

graphlan.py --dpi 150 2.xml 2.png

graphlan.py 2.xml 2.pdf

添加基本全局参数，调图形样式

graphlan_annotate.py --annot annot1.txt 2.xml 3.xml

graphlan.py --dpi 300 3.xml 3.png

graphlan.py 3.xml 3.pdf

问题：分支节点未命名或无法命名

在扩增子分析中，我们经常构建的无根进化树，有没有根倒是无所谓，但是这些使用Fastree构建的进化树有共同的特征，就是除了叶节点，其余分支节点命名无法在graphlan指明。

进化树文件局部示例：

((((((((('OTU_512':0.01655,'OTU_592':0.05858)==1.000:0.15061==,((('OTU_409':0.05856,

但是这里我们可以看到分支节点名为1.000，但是在注释文件中通过这个节点无法添加相应的操作。

参考文献

Asnicar F, Weingart G, Tickle T L, et al. Compact graphical representation of phylogenetic data and metadata with GraPhlAn[J]. PeerJ, 2015, 3: e1029. http://huttenhower.sph.harvard.edu/GraPhlAn

作者简介

最后编辑：2022-06-30

作者：萌小白

一个热爱网络的青年！

站内专栏 邮箱