R: ggtree (二) gheatmap

導(dǎo)讀

上一篇：R語言ggtree繪制進化樹（一）
上一篇相關(guān)：Phylophlan（二）種間進化分析
tree文件和繪圖mapping文件來源：Phylophlan（三）將新物種插入進化樹

一、ggtree畫樹狀圖：長方型、加對齊線

關(guān)鍵參數(shù)：
%<+% 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()

打開繪圖結(jié)果tree_rectangular_line.pdf，如下：

圖片.png

二、ggtree畫樹狀圖：長方型、末端著色

關(guān)鍵參數(shù)：
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()

打開繪圖結(jié)果tree_rectangular_point.pdf，如下：

圖片.png

三、ggtree畫樹狀圖：圓型、加對齊線

關(guān)鍵參數(shù)：
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()

打開繪圖結(jié)果tree_circular_line.pdf，如下：

圖片.png

四、ggtree畫樹狀圖：長方型、末端著色

關(guān)鍵參數(shù);
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()

打開繪圖結(jié)果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 讀樹、豐度讀取合并組

# 準備豐度數(shù)據(jù)
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))
{
    # 每組求和后取平均數(shù)，放在最后一列
    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)

# 樹的切割
# http://www.itdecent.cn/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函數(shù)
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)

旋轉(zhuǎn)

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

gheatmap參數(shù)：

# gheatmap參數(shù)：
offset = 0.6  # heat to tree 偏移
width = 0.5  # heat 寬度
low  # 低值顏色
high  # 高知顏色
color  # cell border 顏色
colnames_level  # 列名label
colnames_angle = 90  # 列名旋轉(zhuǎn)度數(shù)
font.size = 0.1  # colname 大小
hjust=.90  # 列名位置調(diào)整

參考：
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: 系統(tǒng)發(fā)育樹（phylogenetic tree）可視化
7 plotTree-ggtree

\color{green}{????原創(chuàng)文章，碼字不易，轉(zhuǎn)載請注明出處????}

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

  • 導(dǎo)讀 Prodigal是原核生物基因預(yù)測軟件，常被用于原核生物de novo組裝分析中。Prodigal能預(yù)測由d...
    胡童遠閱讀 5,585評論 10贊 15
  • 進化樹的可視化軟件非常多，其中R包 ggtree 功能非常強大，非常靈活，簡單記錄自己的學習筆記 第一步：使用 m...
    小明的數(shù)據(jù)分析筆記本閱讀 16,239評論 12贊 18
  • TaoYan 簡介 本文將繪制靜態(tài)與交互式熱圖，需要使用到以下R包和函數(shù)：heatmap():用于繪制簡單熱圖的函...
    taoyan閱讀 48,019評論 4贊 129
  • 導(dǎo)讀 ggtree由R語言大神Y叔紙筆，于2018年發(fā)表在Molecular Biology and Evolut...
    胡童遠閱讀 26,757評論 2贊 24
  • 圖形初步 在本章中，我們將討論處理圖形的一般方法。我們首先探討如何創(chuàng)建和保存圖形，然后關(guān)注如何修改那些存在于所有圖...
    jplee閱讀 5,325評論 0贊 12

友情鏈接更多精彩內(nèi)容