R: ggtree (二) gheatmap

导读

上一篇：R语言ggtree绘制进化树（一）
上一篇相关：Phylophlan（二）种间进化分析
tree文件和绘图mapping文件来源：Phylophlan（三）将新物种插入进化树

一、ggtree画树状图：长方型、加对齐线

关键参数：
%<+% map # 引入注释文件
layout="rectangular" # 长方型树状图
align=TRUE # 添加对齐虚线
col=Phylum # 以Phylum给虚线着色
legend.position="bottom" # legend至于底部
legend.box="horizontal" # legend水平放置

library(ggplot2)
library(ggtree)
tree=read.tree("all.tree.int.nwk")
data=fortify(tree)
map=read.table("mapping.txt", header=T, sep="\t")

# 长方形（线型）
gra=ggtree(tree, layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=-0.5, align=TRUE, linesize=0.1)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)

pdf("tree_rectangular_line.pdf")
gra
dev.off()

打开绘图结果tree_rectangular_line.pdf，如下：

图片.png

二、ggtree画树状图：长方型、末端着色

关键参数：
ggtree(aes(col=Phylum)) # 末端“枝”颜色
geom_tippoint(aes(color=Phylum)) # 末端“点”颜色

gra=ggtree(tree, aes(col=Phylum), layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端点颜色、大小
geom_tiplab(aes(label=NA), size=0.8)+
# 注释
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)

pdf("tree_rectangular_point.pdf")
gra
dev.off()

打开绘图结果tree_rectangular_point.pdf，如下：

图片.png

三、ggtree画树状图：圆型、加对齐线

关键参数：
layout="circular" # 画圆型树状图

# 圆形（线型）
gra=ggtree(tree, layout="circular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=NA, col=Phylum), hjust=2, align=TRUE, linesize=0.1)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)

pdf("tree_circular_line.pdf")
gra
dev.off()

打开绘图结果tree_circular_line.pdf，如下：

图片.png

四、ggtree画树状图：长方型、末端着色

关键参数;
ggtree(aes(col=Phylum)) # 末端“枝”颜色
geom_tippoint(aes(color=Phylum)) # 末端“点”颜色

gra=ggtree(tree, aes(col=Phylum), layout="circular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tippoint(aes(color=Phylum), size=0.1)+
# 端点颜色、大小
geom_tiplab(aes(label=NA, col=NA), size=0.8)+
# 注释、注释的颜色
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5)))+
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)

pdf("tree_circular_point.pdf")
gra
dev.off()

打开绘图结果tree_circular_point.pdf，如下：

图片.png

五、ggtree画树状图：长方型、加对齐线、加注释

1 准备注释文件

# 制作注释文件
cd /home/cheng/huty/softwares/phylophlan/data
cat ppafull.tax.txt | sed 's/.p__/\t/g' | sed 's/.c__/\t/g' | sed 's/.g__/\t/g' | awk 'BEGIN{OFS="\t"}{print $1, $3, $5}' > tax.txt

# 绘图准备
touch tax_bin.txt
echo -e 'id\tPhylum\tSpecies' > tax_bin.txt

cat bin_taxonomy.txt | sed '1d' | sed 's/.p__/\t/' | sed 's/.c__/\t/' | awk '{printf("%s\t%s\t%s\n", $1, $3, $1)}' >> tax_bin.txt
cat /home/cheng/huty/softwares/phylophlan/data/tax.txt >> tax_bin.txt

2 绘图

# 绘制insert tree
library(ggplot2)
library(ggtree)
tree=read.tree("wm.insert.tree.int.nwk")
data=fortify(tree)

map=read.table("tax_bin.txt", header=T, sep="\t")


gra=ggtree(tree, layout="rectangular", size=0.1) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=Species, col=Phylum), align=TRUE, size=0.18, linesize=0.05)+
# 注释、颜色、高度、对其、虚点大小
theme(legend.title=element_text(face="bold"), legend.position="bottom", legend.box="horizontal", legend.text=element_text(size=rel(0.5))) +
# 图例位置、文字大小
xlim(NA, max(data$x)*1.3)

ggsave(gra, filename="tree_rectangular.pdf", height=48, width=10)

六、ggtree + gheatmap组合图

library(ggplot2) # 加载ggplot2
library(ggtree) # 加载ggtree

tree=read.tree(list.files()[grepl("denovo.tree.nwk", list.files())]) # 读取nwk文件
map=read.table("map.txt", header=T, sep="\t", comment.char="")
data=fortify(tree)

abun = read.table("/home/cheng/client2/wangmeng/X101SC19061144-Z01-J041/Bin_all/Bin_quant/bin_abundance_table.txt", header=T, sep="\t", row.names=1)
abun = abun[, order(colnames(abun), decreasing=F)]

p = ggtree(tree, layout="rectangular", ladderize=FALSE, branch.length="none", size=0.8, aes(color = Phylum)) %<+% map +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(col=Phylum), align=TRUE, size=2.5) +
# 注释、颜色、高度、对其、虚点大小
theme(legend.position = 'none')
ggsave(p, file="bin_pheatmap.pdf", width=8, height=8)

p2 = gheatmap(p, abun, offset = 0.6, width = 0.5, 
             font.size=3, low = "green", high = "red", color = "white", 
             colnames_level = colnames(abun), 
             colnames_angle=90, hjust=.90)

ggsave(p2, file="bin_pheatmap2.pdf", width=14)

用ggtree来绘制进化树

七、ggtree + facet_plot柱形图

1 读树、丰度读取合并组

# 准备丰度数据
abun = read.table("/home/cheng/client2/wangmeng/X101SC19061144-Z01-J041/Bin_all/Bin_quant/bin_abundance_table.txt", header=T, sep="\t", row.names=1)
abun = abun[, order(colnames(abun), decreasing=F)]

# 合并列，样品组
df = abun[, c(1:7)]
k = 1
for(i in c(1, 4, 7, 10, 13, 16, 19))
{
    # 每组求和后取平均数，放在最后一列
    df[, k] = apply(abun[, c(i, i+1, i+2)], 1, sum)/3
    # 最后一列重命名，用组名
    colnames(df)[k] = substr(colnames(abun)[i], 1, 2)
    k = k + 1
}

# tree
tree=read.tree("/media/cheng/disk3/Projects/wangmeng/X101SC19061144-Z01-J041/Bin_all/Bin_phylo/wm.insert/wm.insert.tree.int.nwk")
data=fortify(tree)

# 注释文件
anno = read.table("/media/cheng/disk3/Projects/wangmeng/X101SC19061144-Z01-J041/Bin_all/Bin_phylo/wm.insert/tax_bin.txt", header=1, sep="\t")

2 切割丰度表、切割树

# 制作丰度表
bin = read.table("list2.txt")  ## 挑选bin
table = read.table("table2.txt", header=T, sep="\t", row.names=1)
## 挑选物种
## 0填充丰度

df2=df[rownames(df)%in%bin$V1,] # 提取bin丰度
df2=rbind(df2, table) # 列合并

df2$tax=rownames(df2) # 行名加到注释列
rownames(df2)=1:nrow(df2) # 重排行名

# 注释
anno2=anno[anno$Species%in%df2$tax,]
rownames(anno2) = 1:nrow(anno2)

# 切割丰度表：注释 + 丰度
anno_abun=merge(anno2, df2, by.x='Species', by.y='tax', all.x=T)
ck = data.frame(id=anno_abun$id, num=anno_abun$CK*10)
t1 = data.frame(id=anno_abun$id, num=anno_abun$T1*10)
t2 = data.frame(id=anno_abun$id, num=anno_abun$T2*10)
t3 = data.frame(id=anno_abun$id, num=anno_abun$T3*10)
t4 = data.frame(id=anno_abun$id, num=anno_abun$T4*10)
t5 = data.frame(id=anno_abun$id, num=anno_abun$T5*10)
t6 = data.frame(id=anno_abun$id, num=anno_abun$T6*10)

# 树的切割
# https://www.jianshu.com/p/e1d0a0fce5a2
# ggtree:::getSubtree(tree, data_sub$node)
sub = anno2$id
data_drop = data[!(data$label%in%sub),]
tree_sub = drop.tip(tree, data_drop$label)
data_sub=fortify(tree_sub)

# id 列放在首列
gra = ggtree(tree_sub, aes(color = Phylum),layout="rectangular", size=0.5) %<+% anno2 +
# 树型、线粗细、末端颜色 + 注释信息
geom_tiplab(aes(label=Species, color = Phylum), align=TRUE, size=4, linesize=1) +
# 注释、颜色、高度、对其、虚点大小
theme(legend.position = 'none') +
xlim(0, max(data$x)*2)

ggsave(gra, filename="task2_tree.pdf", height=10, width=10)

3 用facet_plot添加多个柱形图

library("ggstance")
# 需要ggstance提供barh函数
p <- facet_plot(gra, panel = 'CK', data = ck , geom = geom_barh, aes(x=num), fill="#db616a", color="#db616a", stat='identity') %>% 
facet_plot(panel = 'T1', data = t1 , geom = geom_barh, aes(x = num), fill="#dbac61", color="#dbac61", stat='identity') %>%
facet_plot(panel = 'T2', data = t2 , geom = geom_barh, aes(x = num), fill="#6cdb61", color="#6cdb61", stat='identity') %>%
facet_plot(panel = 'T3', data = t3 , geom = geom_barh, aes(x = num), fill="#36bfa2", color="#36bfa2", stat='identity') %>%
facet_plot(panel = 'T4', data = t4 , geom = geom_barh, aes(x = num), fill="#00BFFF", color="#00BFFF", stat='identity') %>%
facet_plot(panel = 'T5', data = t5 , geom = geom_barh, aes(x = num), fill="#836FFF", color="#836FFF", stat='identity') %>%
facet_plot(panel = 'T6', data = t6 , geom = geom_barh, aes(x = num), fill="#EE00EE", color="#EE00EE", stat='identity') +
theme_tree2() +
scale_x_continuous(limits=c(0, 11), breaks = c(0, 5, 10)) +
labs(x="Avg(TPM) X 10") +
theme(legend.position="left", text=element_text(family="serif"))

ggsave(p, filename="task2_tree_bar.pdf", height=10, width=50, limitsize = FALSE)

# combine
hmo = map[, 3:ncol(map)]
rownames(hmo) = map$label
hmo_tmp = hmo
for(i in 1:nrow(hmo))
{
    for(j in 1:ncol(hmo))
    {
        if(hmo[i, j] != 0){hmo_tmp[i, j] = 1}  # 0/1
    }
}

tmp = hmo_tmp
num = seq(10, 100, by=10)
for(i in 1:10)
{
    tmp[,i] = as.character(tmp[,i]*num[i]) # (10, 100, by=10)
}

# base
gra=ggtree(data, color="black", layout="fan", size=1, open.angle=10) %<+% map +
geom_tippoint(aes(col=Phylum), size=5) +
theme(legend.title=element_text(face="bold", size=15), legend.position="right", legend.text=element_text(size=13)) +
theme(text=element_text(family="serif")) +
scale_color_manual(
    values = c("Actinobacteriota" = "#33a02c",
    "Bacteroidota" = "#EE1289",
    "Desulfobacterota" = "#b2df8a",
    "Firmicutes" = "#009ACD",
    "Fusobacteriota" = "#EEAD0E",
    "Proteobacteria" = "#EE6363",
    "Synergistota" = "#EEB4B4",
    "Verrucomicrobiota" = "#9F79EE"))

ggsave(gra, file="tree.pdf", height=30, width=30)

# gheatmap ten color
colours=c("white","#20B2AA","#DC143C","#0000FF","#9370DB","#1E90FF","#7CFC00","#FFFF00","#FF00FF","#FA8072","#7B68EE")
names(colours) <- seq(0, 100, by=10)
cazys = c("GH2","GH20","GH29","GH33","GH42","GH95","GH112","GH136","CBM32","CBM51")

p2 <- gheatmap(gra, tmp, offset = 0, width=0.5, font.size=10, color = NULL, hjust = -0.1, colnames_level=colnames(hmo), colnames_angle=-90, family="serif", legend_title="CAZyme") + scale_fill_manual(values=colours, breaks=seq(0, 100, by=10))

ggsave(p2, file="tree_heat2.pdf", height=30, width=30)

# base
gra=
ggtree(data, color="black", layout="fan", branch.length="none", size=0.5, open.angle=15) %<+% map +
geom_tippoint(aes(col=Phylum), size=5) +
theme(legend.title=element_text(face="bold", size=20),                       legend.position="right",
      legend.text=element_text(size=15),
      legend.key=element_rect(size=20)) +
theme(text=element_text(family="serif")) +
scale_color_manual(
    values = c("Actinobacteriota" = "#33a02c",
    "Bacteroidota" = "#EE1289",
    "Desulfobacterota" = "#b2df8a",
    "Firmicutes" = "#009ACD",
    "Fusobacteriota" = "#EEAD0E",
    "Proteobacteria" = "#EE6363",
    "Synergistota" = "#EEB4B4",
    "Verrucomicrobiota" = "#9F79EE"))

ggsave(gra, file="tree2.pdf", height=30, width=40)

##############################
## gheatmap two color
##############################
colours=c("white","#20B2AA","#20B2AA","#20B2AA","#20B2AA","#20B2AA","#20B2AA","#20B2AA","#20B2AA","#DC143C","#DC143C")
names(colours) <- seq(0, 100, by=10)

##
p3 <- gheatmap(gra, tmp, offset = 0, width=1, font.size=14, color = NULL, hjust = -0.1, colnames_level=colnames(hmo), colnames_angle=-90, family="serif", legend_title="CAZyme") + scale_fill_manual(values=colours, breaks=seq(0, 100, by=10))

ggsave(p3, file="tree_heat2.pdf", height=30, width=30)

旋转

gra2 = rotate_tree(gra, 90)
ggsave(gra2, file="enzyme_tree/tree2.pdf", height=30, width=30)

gheatmap参数：

# gheatmap参数：
offset = 0.6  # heat to tree 偏移
width = 0.5  # heat 宽度
low  # 低值颜色
high  # 高知颜色
color  # cell border 颜色
colnames_level  # 列名label
colnames_angle = 90  # 列名旋转度数
font.size = 0.1  # colname 大小
hjust=.90  # 列名位置调整

参考：
1 ggplot2版聚类物种丰度堆叠图
2 Two Methods for Mapping and Visualizing Associated Data on Phylogeny Using Ggtree
3 ggtree facet_plot gheatmap
4 用ggtree画进化树
5 R doc gheatmap
6 ggtree: 系统发育树（phylogenetic tree）可视化
7 plotTree-ggtree

\color{green}{😀😀原创文章，码字不易，转载请注明出处😀😀}

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 206,968评论 6赞 482
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 88,601评论 2赞 382
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 153,220评论 0赞 344
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 55,416评论 1赞 279
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 64,425评论 5赞 374
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 49,144评论 1赞 285
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,432评论 3赞 401
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,088评论 0赞 261
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 43,586评论 1赞 300
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,028评论 2赞 325
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,137评论 1赞 334
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,783评论 4赞 324
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,343评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,333评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,559评论 1赞 262
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,595评论 2赞 355
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,901评论 2赞 345

推荐阅读更多精彩内容

  • Phylophlan（三）将新物种插入进化树
    导读 Prodigal是原核生物基因预测软件，常被用于原核生物de novo组装分析中。Prodigal能预测由d...
    胡童远阅读 5,132评论 10赞 15
  • ggtree笔记
    进化树的可视化软件非常多，其中R包 ggtree 功能非常强大，非常灵活，简单记录自己的学习笔记 第一步：使用 m...
    小明的数据分析笔记本阅读 15,822评论 12赞 18
  • R语言学习笔记之热图绘制
    TaoYan 简介 本文将绘制静态与交互式热图，需要使用到以下R包和函数：heatmap():用于绘制简单热图的函...
    taoyan阅读 47,390评论 4赞 129
  • R: ggtree (一) base
    导读 ggtree由R语言大神Y叔纸笔，于2018年发表在Molecular Biology and Evolut...
    胡童远阅读 26,154评论 2赞 24
  • R语言图形初步
    图形初步 在本章中，我们将讨论处理图形的一般方法。我们首先探讨如何创建和保存图形，然后关注如何修改那些存在于所有图...
    jplee阅读 4,994评论 0赞 12